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INTRODUCTION 
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y  Optical  surveillance  sensors  are  sometimes  subject  to  static  noise  pro¬ 
cesses  which  can  adversely  affect  the  optimal  processing  of  their  resulting 
imagery.  For  example,  image  intensity  values  derived  from  line  outages,  dead 
pixels,  ‘"popcorn'^noise  and  other  such  noise  mechanisms  will  contaminate  both 
local  and  global  estimates  of  the  power  spectrum  density,  probability  density 
function  and  other  key  statistical  properties  inherent  to  the  original  ob¬ 
served  scene.  If  these  estimates  are  used  to  derive  optimum  filters  for 
detecting  specific  targets  in  said  imagery  or  registering  sequences  of  images, 
poor  processing  performance  could  result.  References  1  through  7  provide 
excellent  reviews  of  current  image  processing  trends  dependent  on  good  quality 
pictures  and  illustrate  their  utility  for  enhancing  the  inherent  information 
content  found  in  remotely  sensed  images  such  as  those  taken  by  the  LANDSAT  and 
NIMBUS-7  satellites.  They  also  show  the  effect  of  noisy  pixels  on  these 
techniques  and  the  types  of  performance  degradations  incurred;  which  can  be 
significant.  This  suggests  that  methods  for  replacing  the  ^’bad"  pixel  values 
with  numbers  commensurate  with  the  inherent  statistics  of  a  detected  image  can 
be  important  to  optimum  image  processing  in  many  applications.  r _ _ 

The  standard  techniques  for  replacing  bad  pixels  is  to  replace  their 
recorded  values  with  the  scene  mean  intensity  value,  or  to  average  the  eight 
intensity  values  surrounding  the  bad  pixel  and  substitute  the  results  for  the 
incorrect  values.  Unfortunately  these  methods  do  not  necessarily  take  into 
account  the  spatial  coherence  properties  of  tne  inherent  surface  clutter  and 
may  not  give  the  best  representation  of  the  true  intensity  values  which  should 
be  there.  The  intent  of  this  report  is  to  describe  the  potential  merits  of 
three  nonlinear  estimation  techniques  for  replacing  noisy  pixels  in  computer 
contaminated  DAEDALUS  imagery.  The  performance  of  these  techniques  will  be 
benchmarked  against  the  neighborhood  averaging  technique  cited  above. 


PROBLEM  APPROACH 


During  the  past  decade  several  nonlinear  estimation  techniques  have  been 
used  to  smooth  spiky  noise-contaminated  time  series  and  optical-image  data, 
and  their  success  has  been  reported  in  references  5,  6,  and  8-10.  Median 
filters  have  emerged  as  one  of  the  best  methods  found  since  they  are  effective 
in  estimating  reasonable  insertion  values,  while  still  preserving  any  mono¬ 
tonic  step  edges  present  in  the  data  (8).  However,  these  filters  and  other 
similar  data  estimators  do  contaminate  the  basic  intensity  statistics  and  this 
can  affect  subsequent  multispectral  or  time  series  data  processing  in  some 
nonlinear  fashion.  This  can  be  especially  detrimental  when  strong  image-to- 
image  correlations  are  desired  (11).  The  question  is,  "which  nonlinear  data 
estimator  provides  the  best  data  representation  after  application?"  In  an 
attempt  to  answer  this  question,  three  popular  nonlinear  data  filtering  tech¬ 
niques  were  assessed  in  their  ability  to  correct  computer  contaminated 
DAEDALUS  imagery  of  significant  outliers,  with  and  without  a  bad  pixel  locali¬ 
zation  routine,  and  create  reasonably  accurate  renditions  of  the  uncontamin¬ 
ated  scenes.  The  benchmark  performance  levels  were  assumed  to  be  those 
obtainable  from  linear  filtering  using  a  neighborhood  average  and  replacement 
technique  on  the  same  noisy  imagery.  In  the  next  two  subsections,  noisy  pixel 
replacement  estimators,  specific  DAEDALUS  imagery,  and  detailed  experimental 
procedures  used  in  this  investigation  will  be  described. 


PROPOSED  NOISY  PIXEL  REPLACEMENT  ESTIMATORS 


Four  data  replacement  estimators  were  assumed  for  this  work;  the  first  a 
linear  technique  and  the  last  three  were  of  a  nonlinear  nature:  namely 

•  Neighborhood  Averaging  Filter 

•  Bi weight  Filter 

•  Median  Filter 

•  Four  way  Median  Filter 


This  subsection  provides  a  detailed  description  of  each.  Each  filter  operates 
on  the  nine  pixels  found  in  a  3x3  pixel  window  about  the  pixel  to  be  replaced 
and  these  pixels  are  denoted  in  the  following  matrix  form: 


x( i— 1, j— 1 )  x  (i- 
x( i, j-1)  x( i , j 
x(i+l,j-l)  x(i+l 


i,j> 

x( 1- 1 , j+1 ) 

) 

x( i , j+ 1 ) 

» j ) 

X (  1+ 1 ,  J  +  1 ) 

In  the  above  matrix,  x(i,j)  is  the  candidate  pixel  to  be  replaced.  For  the 
biweight  and  median  filters,  these  nine  values  were  passed  to  a  sorting 
routine  (SHELLSORT)  so  they  could  be  assembled  in  low  to  high  order.  The 
other  two  filters  did  not  require  sorting  and  no  other  additional  pret'iltering 
operations  were  imposed. 


Neighborhood  Averaging  Replacement  Filter 


The  neighborhood  averaging  replacement  filter  is  one  which  sums  the  ampli¬ 
tude  values  of  the  eight  array  elements  surrounding  x(i,j)  and  divides  the 
result  by  eight  to  yield  the  replacement  value  for  x(i,j).  This  computation 
is  similar  to  a  local  neighborhood  mean  estimate,  except  that  tne  center  ele¬ 
ment  is  missing  from  the  calculation.  The  effectiveness  of  this  estimation 
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method  is  dependent  upon  the  percentage  of  noise  inherent  in  the  image,  the 
spatial  variability  of  the  scene,  and  whether  some  sort  of  noisy  pixel  screen¬ 
ing  procedure  is  applied  to  determine  the  need  for  filtering.  Since  this  is 
the  benchmark  technique  for  the  study,  these  points  will  be  discussed  in 
greater  detail. 

Most  infrared  and  optical  imagery  are  Markov-like  random  processes  and 
have  autocovariance  functions  with  envelopes  which  fall  off  exponentially  over 
some  linear  distance(s).  If  the  clutter  level  of  the  image  is  low,  a  bad 
pixel  replacement  can  result  because  of  either  the  presence  of  outliers,  or 
rapid  decay  of  the  spatial  covariance.  This  latter  situation  implies  the 
spatial  variability  of  the  scene  is  high,  such  as  an  image  of  a  city  or  other 
such  cluttered  terrain.  A  noisy  pixel  test  can  minimize  the  occurrence  of  the 
former,  but  nothing  can  really  be  done  about  the  former  without  more  informa¬ 
tion.  If  the  autocovariance  of  images  falls  off  slowly,  which  implies  slowly 
varying  clutter  fields,  bad  pixel  replacement  will  most  likely  only  occur  with 
outliers  present.  Here  again  a  noisy  pixel  test  will  help  minimize  efforts  in 
pixel  replacement. 

When  the  noise  level  is  high,  bad  pixel  replacements  can  occur  because 
the  mean  and  standard  deviations  derived  from  the  imagery  may  replace  good 
values  with  bad  or  pass  noisy  pixels  off  as  good.  This  is  especially  possible 
when  the  spatial  variability  of  the  clutter  is  high.  Obviously,  a  noisy  pixel 
test  may  not  improve  anything  at  all  if  its  value  cannot  be  easily  determined 
by  the  chosen  test. 

Given  the  above,  nonlinear  estimation  techniques  suffer  the  same  problems 
and/or  to  the  same  degree.  This  is  the  objective  of  the  study  and  will  be  the 
aspect  investigated  in  some  detail  in  the  sections  to  come. 

Biweight  Filter 

For  every  pixel  in  an  image  of  interest,  an  ordered  nine  element  string 
ot  in  tens i t ies , 

x(1 )<x(2)<x(3)...<x(d)<x(9), 

is  [vis  si'ii  tj  t  fie  hi  weight  filter  subroutine  (  BIWEIGHT)  tor  pixel  replacement. 

The  first  operation  performed  is  to  calculate  the  string  standard  deviation 
based  upon  the  interquartile  values  x(3)  and  x(7),  i.e.  the  subroutine 
compu tes 

standard  deviation  (.SO)  =  l  x(  7  )-x(  3  )  ] /I  .  349;  for  x(7).ne.x(3)  (la) 

1 .483  (=1/.674b)  tor  x(7) .eq.x( 3) .  (1b) 


The  subroutine  then  derives  an  estimated  mean  using  those  elements  in  the 
string  within  five  standard  deviations  ot  the  median  element  x(5).  Mathe- 
mifio.iily  t  tie  program  calculates 


23  x(i)*rectlx(i)-x(5))/(10*SD)] 


rect[ (x) ( i )-x( 5 ) )/10*SD) ] 

1 

where  x  is  the  estimated  mean  and  rect(z)  is  the  rectangular  function  given  by 

rect(z)  =  1;  for  abs(z)<1/2 
=  0;  otherwise. 

This  particular  form  of  mean  estimate  allows  one  to  derive  an  estimated  mean 
devoid  of  significant  outliers  which  would  contaminate  the  estimate. 

A  weighting  function  is  then  calculated  using  the  formula 

w(i)=[1-[(x(i)-x)/(6*SD)]**2]**2  (3) 

for  [  (x(  i  )-x) /(6*SD)  ]  **2  >  1.0  the  term  becomes  1.0  and  w(i)  =0.  A  new 

estimate  of  the  mean  is  computed  through  the  relation 


Y  x(i)*w(i) 
i 

9 

Y  w<  i ) 


and  this  number  replaces  the  array  element  x(i,j).  This  technique  has  been 
suggested  as  a  more  robust  means  of  smoothing  noisy  data  sets  than  the  median 
filter  [12]  and  this  assertion  will  be  evaluated  in  the  study. 

Median  Filter 

For  every  pixel  in  an  image  of  interest,  an  ordered  nine  element  string 
of  intensities, 


x(1 )<x(2)<x(3)...<x(8)<x(9), 

is  passed  to  the  median  filter  subroutine  and  the  median  string  element  x(5) 
is  substituted  for  x(i,j).  The  major  advantage  of  median  filtering  is  that 
constant  backgrounds,  slopes  and  edges  are  preserved,  while  isolated  pulses 
less  than  or  equal  to  m=(n-1)/2  are  suppressed,  n  being  the  lenqth  of  the  data 
string  and  equal  to  9  in  this  case. 

Four-way  Median  Filter 

The  four-way  median  filter  is  a  filtering  routine  which  applies  a  one¬ 
dimensional  median  filter  in  the  north-south,  east-west,  southwest-northeast 
and  northwest-southeast  directions  sequentially.  Pictorial ly,  the  filter 
operates  in  the  following  sequence  order: 


Although  this  filter  only  operates  on  three-pixel  strings  at  any  one  pass, 
this  pixel  estimation  technique  requires  a  larger  initial  array  size  upon 
which  to  operate  than  the  previous  techniques.  This  is  because  the  subroutine 
performs  an  in-place  four-way  filtering  operation  and  the  resulting  array  size 
from  each  pass  is  reduced  by  2  rows  and  2  columns. 

Unlike  the  neighborhood  average  and  biweight  filters  previously  de¬ 
scribed,  this  method  does  not  replace  array  elements  with  computed  values,  but 
rather  replaces  them  with  one  of  the  nine  original  pixel  values.  The  replace¬ 
ment  may  be  a  good  or  noisy  pixel,  depending  on  its  relationship  to  its  adja¬ 
cent  neighbors  on  each  pass,  and  whether  a  bad  pixel  test  is  applied.  Tech¬ 
nically,  this  type  of  filter  derives  a  proper  replacement  only  when  one  noisy 
pixel  exists  in  any  chosen  3x3  pixel  array  and  may  make  an  error  otherwise  if 
more  than  one  outlier  exists  in  these  windows.  Clearly  the  particular  posi¬ 
tion  of  the  noisy  pixels  relative  to  the  order  sequence  of  the  filter  will 
define  whether  an  error  is  made  or  not.  For  example,  all  the  pixels  in  the 
3x3  matrix  could  be  noisy  except  for  the  upper  left  and  tne  lower  right  cor¬ 
ners.  On  the  fourth  pass,  a  good  replacement  would  occur. 

TEST  IMAGERY 

The  test  imagery  chosen  for  this  investigation  is  a  set  of  512x512x8  bits 
DAEDALUS  derived  images  provided  by  E.M.  Winter  of  Technical  Research  Asso¬ 
ciates  [13].  The  DAEDALUS  sensor  was  flown  in  the  NASA/AMES  U-2  over  San 
Francisco  Bay  and  points  south  on  14  September  1983,  and  specific  locations 
interrogated  included  San  Jose,  Santa  Cruz,  Monterey  Bay  and  Point  Mugu, 
California.  The  images  used  in  this  study  were  from  the  channel  5  segment  and 
possess  an  individual  pixel  field-of-vzew  of  approximately  28  meters.  More 
information  on  this  flight  can  be  obtained  by  contacting  Dr.  Winter  at  (805) 
987-1972  in  care  ot  Technical  Research  Associates,  Inc.,  445  Rosewood  Avenue, 
Suite  H,  Camarillo,  California  93010. 

'this  set  was  chosen  because  of  the  diversity  of  terrain  recorded,  i.e., 
the  range  of  spatial  variability  interrogated.  The  six  images  selected  tor 
the  investigation  were: 

1.  Lower  dan  Francisco  Bay  (SKI):85%  bay,  15%  city  and  tlatlands 

2.  Northern  dan  Jose  (SK2):85%  city,  15%  bay 

3.  Downtown  dan  Jose  (SK3):95%  city,  5%  tlatlands 

4.  Santa  Cruz  Mountains  (SF6):10U%  mountains 

5.  Santa  Cruz  (SF8):30%  city,  7u%  ocean 

h  • 


ocean  and  clouds  oft  Santa  Cruz  (SF11):30%  ocean,  70%  clouds 
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These  images  are  shown  in  their  uncontaminated  form  in  appendix  B.  In  the 
work  to  be  described,  two  versions  of  tne  above  imagery  were  used.  One  ver¬ 
sion  used  the  calibrated  DAEDALUS  imagery  directly  in  the  study.  Unfortu¬ 
nately,  most  of  the  images  possessed  pixel  values  between  30  and  120,  and  tnis 
gave  a  very  low  contrast  image  on  the  AED512  displays  used.  To  alleviate  this 
problem,  a  second  version  was  created  which  had  the  six  images  minimum/maximum 
scaled  into  an  8  bit-digitized  intensity  form,  i.e.,  the  minimum  value  of  the 
image  was  set  to  0,  and  the  maximum  value  set  to  255  and  the  rest  of  the  pixel 
values  were  sorted  into  the  254-integer  levels  in  between.  This  allowed  bet¬ 
ter  definition  on  the  display  and  a  better  sense  of  reality  to  tne  viewer. 
However,  it  does  affect  the  inherent  image  statistics  in  a  nonlinear  way. 

Figure  1  illustrates  the  resulting  intensity  histograms  for  lowe'-  San 
Francisco  Bay  scene.  It  is  apparent  that  the  figure  possesses  very  similar 
envelopes,  as  one  would  expect.  The  scaled  DAEDALUS  histogram  is  simply  an 
expanded  version  of  the  unsealed  DAEDALUS  by  a  factor  2.83.  All  of  tne 
DAEDALUS  imagery  used  in  this  study  exhibited  similar  type  scaling.  The  im¬ 
pact  of  this  scaling  is  that  the  scaled  images  will  produce  much  larger  differ¬ 
ence  variances  than  the  unsealed  images  when  both  are  subjected  to  the  same 
statistical  manipulation.  However,  this  is  an  artificial  difference.  The 
relative  performance  of  the  various  filters  on  the  two  image  versions  is  the 
same  and  one  will  be  able  to  draw  the  same  conclusions  by  focusing  upon  one  or 
the  other  set.  In  other  words,  absolute  variances  derived  from  unsealed 
results  should  not  be  compared  with  those  obtained  from  scaled  results. 
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Figure  1.  Intensity  histograms  fot  raw  and  scaled  lower  San  Francisco  Bay 
DAEDALUS  images. 
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EXPERIMENTAL  PROCEDURE 


Programs  (see  appendix  C)  were  written  to  allow  the  operator  to  select  an 
image,  read  it  into  memory,  offset  the  entire  image  from  zero  by  +300,  and 
allow  selection  of  the  following  options: 

V. 

•  Introduction  of  Gaussian  noise:  percent,  mean,  and  sigma  of  noise 

•  Bad  pixel  check  (skip  or  perform) 

•  Number  of  sigmas  from  the  mean  of  the  selected  neighborhood  (for  bad 
pixel  check) 

•  Neighborhood  size;  3x3,  5x5,  7x7,  9x9,  11x11  (for  bad  pixel  check) 

After  selection  of  these  options,  the  program  introduced  noise  or  not, 
then  displayed  the  mean  and  sigma  of  the  difference  between  the  original  and 

noise  contaminated  images.  The  program  then  either  skipped  or  performed  the 

bad  pixel  check  on  a  452x452  subset  of  the  original  image.  If  the  bad  pixel 
check  was  desired,  the  mean  and  sigma  of  the  selected  neighborhood  around  the 
pixel  was  compared  with  the  value  of  the  pixel  itself,  and  if  it  was  within 

the  selected  number  of  sigmas  from  the  neighborhood  mean,  it  was  determined  to 

be  a  good  pixel,  and  would  not  have  to  be  filtered. 

If  the  bad  pixel  check  was  skipped,  then  every  pixel  was  filtered;  when 
the  bad  pixel  check  was  performed,  only  the  outlier  pixels  had  the  filtering 
technique  applied. 

After  the  application  of  the  selected  filtering  algorithm,  the  program 
displayed  the  correlation  coefficient  between  the  original  and  the  filtered 
images,  the  mean  and  sigma  of  the  difference  between  the  original  and  the 
filtered  images,  the  percent  of  noise  contamination  selected,  and  the  mean  and 
sigma  of  the  added  noise.  In  addition,  the  numbers  of  pixels  replaced  by 
noise  and  the  number  of  pixels  corrected  by  the  filter  were  displayed.  The 
number  of  noise  pixels  corrected  and  the  number  of  good  pixels  disturbed  were 
also  displayed. 

NOISE  REPLACEMENT 

The  procedure  for  noise  replacement  was  first  to  move  the  range  of  the 
pixel  values  from  0-255  to  300-555  for  the  scaled  DAEDALUS  images,  and  from 
approximately  30-120  to  330-420  for  the  original  DAEDALUS  images.  This  was 
done  by  adding  300  to  every  pixel  value  after  the  image  was  taken  from  disk. 
By  this  method,  the  introduction  of  Gaussian  noise  with  a  mean  of  0.0  and  a 
sLgma  of  5.0  would  make  the  noise  distinctive  from  the  original  pixel  values. 
Noise  values  replacing  original  image  values  were  then  -20  to  +20,  with  the 
majority  close  to  zero. 


7 


The  actual  noise  replacement  was  done  randomly  over  a  462x462  subset  of 
the  512x512  image.  This  allowed  the  actual  experimental  window  of  452x452, 
over  which  the  bad  pixel  check  and  the  filtering  would  occur,  to  be  well  with¬ 
in  the  noise-contaminated  area.  For  the  four-way  median  filter,  the  contami¬ 
nated  area  was  slightly  larger,  468x468,  due  to  the  added  size  of  the  image 
required,  as  the  filter  shrinks  the  image  by  two  rows  and  columns  for  each  of 
the  four  filter  passes,  resulting  finally  in  the  standard  452x452  image.  As  a 
result,  the  percentage  of  noise  replacement  appears  to  be  less  than  that  of 
the  other  filtering  algorithms;  on  the  other  hand  since  the  four-way  median 
filter  makes  four  passes;  the  percentage  of  noise  dealt  with  is  equal  in  per¬ 
cent  to  the  noise  dealt  with  by  the  other  filtering  methods. 

As  the  program  was  repeated  for  different  percentages  of  noise  and  with 
and  without  the  bad  pixel  check,  and  for  six  different  images,  the  random 
stream  was  always  the  same  so  that  no  one  filter  was  given  a  different,  or 
perhaps  adverse  set  of  noise  upon  which  to  work. 

The  program  allows  replacement  with  Gaussian  means  and  sigmas  other  tnan 
0.0  and  5.0.  The  decision  to  use  these  values  was  based  upon  the  need  to  use 
noise  that  was  different  from  the  original  image  pixel  values,  and  thus  stand 
out. 

Noise  replacement  values  used  in  the  experiment  were  0,  1,  2,  3,  4,  5,  6, 
7,  8,  9,  10,  20,  and  30  percent.  Zero  percent  was  selected  so  that  a  gauge 

could  be  established  on  how  much  the  particular  filter  would  disturb  the  image 
when  there  was  no  noise  replacement.  By  careful  examination  of  the  original 
DAEDALUS  images,  flaws  were  discovered  in  the  scanning  sensor  output  in  every 
image.  Therefore,  it  is  reasonable  to  assume  that  in  every  filtering  action, 
some  of  this  "noise"  would  be  replaced. 

BAD  PIXEL  CHECK 


The  program  allows  for  two  variations  in  the  bad  pixel  check;  first,  to 
skip  it  completely,  and  second  to  vary  the  number  of  sigmas. 


If  the  bad  pixel  check  is  skipped,  then  the  filter  is  applied  to  every 
pixel  of  the  image,  which  for  a  452x452  image  is  204,304  pixels.  This  is  the 
maximum  work  for  the  filter,  and  the  test  of  how  well  or  badly  ttie  filter 
performs  in  correcting  the  noisy  pixels  and  not  disturbing  too  many  of  the 
good  pixels. 


The  second  variation,  the  number  of  sigmas,  determines  how  much  work  the 
filter  does.  During  the  bad  pixel  check,  the  sigma  and  mean  are  computed  for 
the  selected  neighborhood  (3x3,  5x5,  7x7,  9x9,  or  11x11)  around  tne  pixel 

which  is  a  candidate  for  filtering.  If  the  candidate  pixel  is  within  plus  or 
minus  the  selected  number  of  sigmas  from  the  mean,  the  pixel  is  considered 
good  and  does  not  require  filtering.  If  it  falls  beyond  the  selected  number 
of  sigmas  about  the  mean,  then  the  pixel  is  a  candidate  for  filtering.  By 
reducing  the  number  of  sigmas,  more  filtering  is  achieved;  by  increasing  the 
number  of  sigmas;  filtering  is  reduced. 


If  the  filtering  is  reduced,  less  harm  is  done  when  the  candidate  pixel 
is  one  of  the  original  image  pixels,  but  when  the  candidate  pixel  is  a  noisy 
pixel,  then  it  may  escape  detection  as  a  bad  pixel  when  the  number  of  sigmas 
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is  high  (i.e.,  3.0  or  more).  On  the  other  hand,  when  the  number  of  sigmas  is 
low  (making  the  good  pixel  window  smaller)  more  pixels  will  be  filtered;  some 
good  and  some  of  the  noisy  pixels.  This  is  the  tradeoff  in  using  the  bad 
pixel  check;  a  number  of  sigmas  has  to  be  selected  which  will  detect  the  noisy 
yet  not  filter  the  good  pixels.  Several  of  the  filters,  by  nature  of  their 
design,  tend  to  do  less  harm  to  good  pixels,  and  have  more  success  with  elimi¬ 
nation  of  noisy  pixels.  These  are  the  median  and  four-way  median  filters. 
These  filters  do  not  arithmetically  modify  pixel  values;  but  by  the  process  of 
elimination  move  pixels  around.  Noise  pixels  become  the  outliers  and  pixels 

within  the  neighborhood  of  the  candidate  pixel  take  the  place  of  the  noisy 

pixel.  The  neighboring  pixels  are  not  the  same  as  the  replaced  pixel,  but  are 
very  close  to  its  value  if  the  gradient  of  the  neighborhood  is  low. 

BAD  PIXEL  CHECK  NEIGHBORHOOD 

Another  variable  tor  the  bad  pixel  check  is  the  neighborhood  of  the  candi¬ 
date  pixel.  The  operator  is  given  the  opportunity  of  selecting  3x3,  5x5,  7x7, 
9x9,  or  llxll  submatrix  for  the  neighborhood.  During  the  bad  pixel  check,  the 
pixels  in  the  submatrix  selected  are  averaged,  while  at  the  same  time  the 
standard  deviation  of  the  neighborhood  is  computed.  Then  the  candidate  pixel 
is  checked  to  see  if  it  is  within  plus  or  minus  the  number  of  sigmas  around 

the  mean.  It  so,  the  pixel  is  considered  good;  if  not,  it  is  filtered.  The 

size  of  the  neighborhood  affects  the  outcome  of  the  bad  pixel  check;  if  the 
neighborhood  is  small,  one  or  two  bad  pixels  in  the  neighborhood  can  contrib¬ 
ute  to  the  mean  and  standard  deviation  in  such  a  way  as  to  bias  the  true  mean 
of  the  neighborhood .  This  is  shown  in  figure  2. 


Figure  2.  Effect  of  noise  on  image  statistics. 


For  the  zero-noise  case,  the  mean  and  sigma  for  the  3x3  neighborhood  are 
430  and  22;  after  introduction  of  one  noisy  pixel,  the  mean  dropped  to  383, 
and  the  sigma  increased  to  137.  On  the  other  hand,  for  the  11x11  neighbor¬ 
hood,  the  mean  shifted  down  from  492  to  488,  and  the  sigma  increased  from  120 
to  127,  with  the  introduction  of  only  one  noisy  pixel.  This  information  was 
produced  with  simulated  data  for  an  image  which  had  a  moderate  gradient.  The 
main  purpose  of  including  the  information  is  to  show  the  effects  of  small 
amounts  of  noise  on  small  neighborhoods,  and  indicate  how  a  single  bad/noisy 
pixel  can  radically  affect  the  normal  statistical  measures  of  quality  on  appli¬ 
cations  in  image  processing. 

In  figure  3,  the  combined  effects  of  number  of  sigmas  and  neighborhood 
size  selections  are  shown  upon  the  numbers  of  noisy  pixels  and  good  pixels 
replaced.  This  study  was  made  with  only  one  percent  noise  using  the  median 
filter  on  the  ocean  and  cloud  image  (SF11DED).  The  plots  3x3g  and  3x3b  repre¬ 
sent  respectively  the  amount  of  good  and  had  pixels  replaced  by  the  action  of 
the  filter  with  increasing  number  of  sigmas  selected  for  the  3x3  neighborhood. 
The  3x3b  plot  at  number  of  sigmas  =  2  replaces  less  than  the  total  (2076) 
noisy  pixels,  so  its  effectiveness  is  not  perfect  at  number  of  sigmas  =2.  In 
addition,  at  number  of  sigmas  =  2,  the  3x3  (3x3g)  neighborhood  is  responsible 
for  replacing  (and  perhaps  disturbing)  4000  good  pixels.  At  number  of  sigmas 
=  1  (not  shown),  the  3x3  neighborhood  replaces  all  of  the  2076  noisy  pixels, 
but  disturbs  more  than  61  ,000  good  pixels.  The  best  operating  point  for  the 
3x3  neighborhood  is  at  number  of  sigmas  =  2.8265,  where  1906  bad  pixels  are 
replaced  and  only  104  good  pixels  are  disturbed.  The  neighborhood  which  per¬ 
forms  the  bad  pixels  check  best  is  the  9x9  neighborhood  with  number  of  sigmas 
=  4,  where  2051  noisy  pixels  are  replaced  and  only  six  good  pixels  are  dis¬ 
turbed.  Of  course  these  plots  were  made  for  the  most  benign  (least  cluttered) 
image  in  the  study  ( SF1 1 DED . IMG ) . 
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The  apparent  result  of  this  analysis  is  to  select  a  large  neighborhood. 
The  original  set  of  data  (shown  in  appendix  A)  was  taken  with  the  bad  pixel 
check  using  a  5x5  neighborhood  and  number  of  sigmas  =  2.  When  the  discovery 
was  made  that  using  the  bad  pixel  check  with  the  5x5  neighborhood  and  the 
number  of  sigmas  =  2  resulted  in  not  detecting  all  of  the  bad  pixels,  and  at 
the  same  time  disturbing  many  good  pixels,  more  data  were  taken  using  an  11x11 
neighborhood  and  varying  numbers  of  sigmas  depending  upon  the  amount  of  noise 
in  the  image.  The  number  of  sigmas  varied  from  4.1  for  zero  noise  to  1.0  for 
40%  noise. 


SIMULATION  GROUND  RULES  AND  LIMITATIONS 


This  section  outlines  the  major  ground  rules  and  limitations  of  the  per¬ 
formance  and  are  described  in  the  following  sections. 

DIGITAL  PROCESSING  PROCEDURES 

In  the  initial  design  of  the  experiment,  everything  was  done  to  assure 
uniformity  of  procedures  for  each  filter,  so  that  comparison  of  filtering 
techniques  would  be  based  upon  identity  of  noise-contaminated  images  for  each 
filter  input.  The  random  number  generator  was  seeded  at  the  same  starting 
point  so  that  every  image  for  each  filter  would  have  the  same  value  of  noise 
entered  into  the  same  pixels.  The  selection  of  a  452x452  submatrix  image  of 
the  original  512x512  DAEDALUS  images  was  standard  throughout  the  study.  The 
introduction  of  noise  was  over  a  462x462  submatrix  for  the  biweight,  median, 
and  neighborhood  replacement  filters.  The  four-way  median  filter  required 
noise  to  be  introduced  over  a  468x463  submatrix,  due  to  the  shrinking  of  the 
four-way  median  filter  through  each  of  its  four  passes,  resulting  finally  in 
the  452x452  image  which  was  compared  to  the  original  image.  The  geometry  for 
the  submatrix  sizes  is  shown  in  figure  4. 


In  the  design  of  the  programs,  options  were  added  to  vary  the  mean  and 
sigma  of  the  Gaussian  noise  to  be  replaced  in  the  image  and  to  vary  the  number 
of  sigmas  for  the  bad  pixel  check.  In  addition,  provision  was  made  tor  the 


selection  of  any  one  of  the  five  different  neighborhoods  for  the  bad  pixel 
check.  These  were  extra  features  which  extend  the  capability  of  testing. 

DESIGN  LIMITATIONS 

Some  of  the  limitations  caused  by  program  design  considerations  are  the 
images  themselves.  Initially  taken  from  the  DAEDALUS  tapes,  these  images  are 
in  raw,  byte  form,  eight  bits  per  pixel.  The  raw  data  image  pixels  have  a 
limited  integer  range,  usually  from  30  to  120.  In  order  to  display  these 

images  and  obtain  a  good  definition  of  geographic  features,  the  programmer 
must  scale  the  images  into  the  color  table  (or  gray-scale)  range  of  0  to  255. 
Both  types  of  images,  the  raw  and  the  scaled  images,  were  included  in  this 
study.  The  difficulty  of  using  various  filters  is  that  floating  point  values 
generated  by  the  filtering  algorithm  (for  example,  the  biweight  and  the  neigh¬ 
borhood  replacement  filters)  are  rounded  to  the  nearest  integer  before  being 

placed  back  into  the  image.  When  the  difference  between  the  original  and  the 

filtered  images  is  taken,  these  roundoff  errors  have  been  insignificant  due  to 
the  randomness  of  the  dropped  fraction  and  the  fact  that  NINT  (nearest  inte¬ 
ger)  function  is  used.  If  the  pixel  value  is  nn.00  to  nn.4999,  only  those 

amounts  are  lost.  If  the  pixel  value  is  nn.5000  to  nn.9999,  the  corresponding 

amounts  are  gained.  Thus,  over  a  large  n,  the  differences  disappear.  This 

has  been  verified  by  a  program.  The  differences  are  not  visible  until  the 
third  decimal  place,  and  the  differences  are  inconsistent  depending  upon  the 
amount  of  noise  induced  into  the  image. 

In  the  illustrations  of  the  filtered  images,  a  small  border  of  noise 

remains  around  the  edge  of  the  image.  This  five-pixel  noise  border  is  outside 

of  the  452x452  filtered  images  area,  but  noise  has  been  added  to  this  area 

(462x462)  to  simulate  a  full  image  with  noise.  The  five-pixel  border  allows 
for  an  11x11  neighborhood  for  the  bad  pixel  check.  The  noise  area  for  the 

tour-way  median  filter  also  contains  a  five-pixel  border  with  residual  noise 
tor  the  same  reason.  The  four-way  median  filter  five-pixel  noise  border  is 
three  pixels  in  each  direction  farther  out  from  the  other  three  filters'  noise 
border . 

FILTER  WEAKNESSES 

In  the  original  study,  the  et tecti veriess  of  the  filters  begin  to  deteri¬ 
orate  rapidly  after  10%  noise.  Beyond  30%,  the  filters  begin  to  lock  onto 

noise  and  begin  to  reject  the  remaining  good  pixels.  Tests  (not  shown)  have 
been  conducted  to  venf  this  condition.  However,  it  is  not  expected  that 

images  with  more  than  10%  noisy  pixels  would  be  worth  processing. 

The  complexity  and  consequent  time  involved  in  running  these  filters  on  a 
computer  are  not  too  great;  but  compared  among  themselves,  there  is  quite  a 
difference  in  execution  time. 

Trie  simplest  is  the  neighborhood  replacement  filter.  It  requires  no 
sorting  or  multiple  passes.  It  simply  sums  the  eight  neighbors  and  divides  by 
eight . 

The  next-  in  complexity  is  the  median  filter,  which  assembles  the  nine 
pixels  from  the  3x3  array  around  the  candidate  pixel  into  a  1x9  array.  This 
array  is  then  passed  to  the  routine  SHELLSORT,  which  on  the  average  makes  six 


passes  to  sort  the  nine  numbers.  The  routine  would  be  more  efficient  if  there 
were  more  than  nine  values.  When  the  array  passed  to  SHELLSORT  is  already  in 
low  to  high  order,  three  passes  are  required  which  perforins  20  IF  statements, 
but  no  movement  of  values.  When  the  routine  is  done,  the  array  contains  the 
values  in  ascending  order.  The  last  step  of  the  filter  is  to  store  the  median 
value  x(5)  into  the  candidate  pixel. 

A  complicated  filter  is  the  biweight  filter.  It  performs  the  same  steps 
as  the  median  filter,  assembling  the  1x9  array,  passing  it  to  SHELLSORT,  but 
on  return  from  SHELLSORT,  the  biweight  routine  is  called  to  do  the  rest  of  the 
processing  as  described  above.  Formerly,  biweight  was  an  iterative  routine 
which  refined  the  answer  through  as  many  as  200  iterations.  But  since  the 

test  for  convergence  was  0.01  between  the  former  (xhat)  and  the  new  (xhat), 

this  portion  of  the  routine  was  removed.  The  operation  of  the  NIN'f  statement 
in  converting  the  floating  point  result  to  integer,  more  than  compensated  for 
any  further  refinement. 

The  application  of  the  four-way  median  filter  requires  four  passes, 
making  it  the  most  time-consuming  filter  in  this  study.  The  actual  method 

used  was  to  do  an  in-place  four-way  filter  reducing  the  loop  control  by  two 
rows  and  two  columns  each  iteration.  Since  tor  each  iteration  the  number  of 
pixels  involved  is  only  tnree,  the  sort  routine  was  not  used.  In  figure  5, 

five  IF  statements  were  used  to  solve  the  truth  table. 

Programming  application  requires  that  when  applying  a  two-dimensional 
loop  control  to  an  nxn  matrix,  the  solution  of  a  matrix  operation  proceeds  row 
by  row  and  column  by  column.  When  in-place  filtering  is  performed  using  a  3x3 
window,  as  all  the  selected  filters  in  this  study  are,  the  following  phenomena 
becomes  apparent  as  snown  in  figure  b. 

The  pixels  above  the  bold  line  are  pixels  which  have  already  been  fil¬ 

tered;  the  pixels  below  the  bold  line  are  pixels  which  have  not  been  filtered, 
including  the  candidate  pixel.  This  is  an  example  of  in-place  filtering. 

There  are  problems  and  also  advantages  to  this  method  or  processing.  The 
advantages  are: 

1.  If  noise  nas  been  eliminated  from  tne  pixels  that  have  already  oeen 
processed,  the  neighbor nood  mean  and  sigma  tor  the  bad  pixel  check  wiLl  aid  in 
determining  more  accurately  whether  the  canlidate  pixel  is  good  or  has  to  be 

f i 1 tered . 

2.  If  the  pixel  has  to  be  filtered,  the  filtering  algorithm  will  nave- 

had  fewer  noisy  pixels  to  deal  with  in  determinin')  a  substitute  value. 

The  problems  are: 

1.  If  a  good  pixel  has  been  replaced  by  another,  or  a  noisy  pixel  has 
not  been  filtered  and  still  remains  m  the  already  processed  area,  then  the 
neighborhood  mean  and  sigma  for  the  bail  pixel  check  will  possibly  result  m 

tagging  a  good  pixel  bad,  or  tagging  a  noisy  pixel  as  good.  In  either  case, 
this  is  not  desirable. 


2.  In  the  filtering  algorithm,  whether  the  candidate  pixel  is  good  or 
bad,  the  chances  are  that  it  wlII  be  replaced  with  one  that  is  either  worse  or 
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almost  as  bad.  Unfortunately,  some  of  the  advantages  and  some  of  the  problems 
prevail  throughout  the  filtering  techniques  studied  in  this  experiment. 

The  solution  then,  could  be  to  not  perform  in-place  filtering.  This  is 
accomplished  by  building  the  filtered  image  into  a  separate  image.  This  was 
tried,  and  in  all  cases  the  results  were  slightly  degraded,  which  means  that 
the  advantages  of  in-place  filtering  slightly  out-weigh  the  disadvantages,  of 
course  when  the  bad  pixel  check  is  skipped,  a  large  part  of  the  problems  of  in- 
place  filtering  disappears,  and  the  advantages  become  preponderant. 

Experimentation  was  performed  with  a  noisy  mean  other  than  0;  a  mean  of 
50  was  tested  which  made  all  of  the  noise  positive,  and  as  expected,  the  dif¬ 
ference  sigmas  were  slightly  improved. 

The  closer  the  mean  of  the  noisy  pixels  approached  the  mean  of  the  good 
image  pixels,  the  difficulty  in  detecting  the  noisy  pixels  become  more  pro¬ 
nounced.  As  shown  in  figure  7,  the  upper  tail  of  the  noisy  pixel  distribution 
intrudes  upon  the  lower  tail  of  the  good  pixel  distribution.  When  this  hap¬ 
pens,  a  lower  number  of  sigmas  for  the  bad  pixel  check  results  in  a  smaller 
number  of  noisy  pixels  replaced  and  a  greater  number  of  good  pixels  disturbed. 
The  number  of  sigmas  for  the  bad  pixel  check  must  include  all  the  good  pixels, 
while  at  the  same  time  excluding  all  of  the  noisy  pixels. 


PIXEL  INTENSITY 


Figure  7.  Effect  of  increasing  noise  mean  on  bad  pixel  detection. 


RESULTS 

In  addition  to  the  filter  performance  results,  this  study  provided  find¬ 
ings  regarding: 

1.  The  use  of  raw  versus  scaled  DAEDALUS  images. 

2.  The  use  of  the  bad  pixel  check  versus  applying  the  filter  over  the 
entire  image. 

3.  Characteristics  of  increasing  noise  and  the  effects  upon  the  filters' 
performances . 

4.  Measurement  statistic  for  filter  performance. 

5.  Filtering  of  different  images. 

6.  Maximum  performance  expected  from  a  given  filter. 

7.  Image  distortion  with  zero-noise. 

8 .  Improvement  of  the  bad  pixel  check. 

RAW  VERSUS  SCALED  DAEDALUS  IMAGES 

In  all  of  the  raw  DAEDALUS  image  plots,  all  four  filters  performed  better 
than  when  the  scaled  DAEDALUS  images  were  used.  The  performance  was  similar 
for  both  types  of  data,  but  the  raw  images  resulted  in  lower  difference  sigmas 
between  the  original  and  filtered  versions.  This  was  due  to  the  narrower 
range  of  histogram  values  exhibited  by  the  raw  DAEDALUS  data  which  were  usual¬ 
ly  between  10  and  120,  as  opposed  to  the  scaled  DAEDALUS  data  which  were  be¬ 
tween  0  and  255.  These  differences  appear  significant,  but  they  are  only 
relative . 

BAD  PIXEL  CHECK  VERSUS  ENTIRE  IMAGE  FILTERING 

The  ranking  of  the  filters  changed  radically  when  the  filter  was  applied 
over  the  entire  image  (no  bad  pixel  check),  and  specifically,  the  neighborhood 
replacement  filter  performance  declined  adversely  over  the  low  range  of  noise 
(0  to  4  percent),  see  figure  8. 

Generally  using  the  bad  pixel  check,  the  following  results  were  true: 

1.  At  0  and  1  percent  noise  the  bi weight  filter  was  better  than  the  rest 
sometimes  by  only  a  very  slight  mar  (in.  The  four-way  filter  outperformed  the 
biweight  it  0  percent  in  3  out  of  the  >3  images,  doing  better  on  land  alone,  as 
opposed  to  the  biweight's  better  performance  when  water  was  involved  in  the 
image . 

2.  At  7  percent  noise,  sometimes  the  biweight,  median,  or  four-way  was 
the  best.  With  the  raw  DAEDALUS  mages,  the  biweight  didn't  do  as  well;  a  tie 
usually  resulted  between  the  median  and  the  four-way  median. 
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3.  At  3  to  5  percent  noise,  the  four-way  filter  was  usually  the  best, 
except  occasionally  the  median  filter  outperformed  the  four-way  filter  using 
the  scaled  versions  of  the  images. 

4.  From  6  to  30  percent  the  four-way  filter  outperformed  all  the  other 
filters  for  both  the  scaled  and  raw  images. 

Generally,  not  using  the  bad  pixel  check  (when  the  filter  was  applied 
over  the  entire  image)  the  following  results  were  true: 

1.  At  0  and  1  percent  noise,  the  biweight  filter  performed  the  best  for 
both  raw  and  scaled  images,  but  only  by  a  slight  amount  over  the  median  filter. 

2.  At  2  percent  noise  the  median  filter  performed  best  for  the  raw 
images,  and  the  biweight  filter  performed  best  for  the  scaled  images.  In  the 
case  of  the  ocean  scene  the  four-way  filter  performed  best  for  the  raw  image. 

3.  At  3  to  30  percent  noise,  the  median  filter  performed  best  for  four 
of  the  scaled  images  (bay,  downtown  San  Jose,  mountains,  and  Santa  Cruz).  The 
four-way  filter  shared  best  performance  with  the  median  filter  in  this  range 
for  two  of  the  scaled  images  (North  San  Jose  and  ocean)  .  For  the  raw  images 
(San  Francisco  Bay  and  downtown  San  Jose),  the  median  filter  performed  best. 
For  four  raw  images  (North  San  Jose,  mountains,  Santa  Cruz,  and  ocean)  the 
four-way  filter  shared  with  the  median  filter,  in  this  range. 

The  results  of  best  filter  performance  are  shown  in  table  I. 

EFFECTS  OF  NOISE  UPON  FILTERING  ACTION 

Figure  9  shows  the  effects  of  increasing  noise  upon  pixel  replacement  for 
the  bad  pixel  check.  Though  this  plot  is  for  the  ocean  image,  the  remaining 
images  have  very  similar  plots.  The  data  for  the  no  bad  pixel  check  (when  the 
filter  is  applied  over  the  entire  image)  are  not  worth  considering  since  the 
number  of  pixels  replaced  are  204,304  for  the  biweight,  median,  and  neighbor¬ 
hood  replacement  filters,  and,  though  the  numbers  vary,  between  431,000  and 
619,000  for  the  four-way  filter. 

The  interesting  aspect  of  figure  9  is  that  it  clearly  depicts  the  actions 
of  the  filters  throughout  the  noise  range  of  0  to  30  percent.  The  key  points 
are  0  and  10  percent.  At  0  percent  noise,  all  of  the  filters  perform  filter¬ 
ing  when  in  reality,  none  is  required.  This  of  course  means  that  this  ten¬ 
dency  of  a  filter  to  replace  good  pixels  is  going  to  persist  to  some  degree  as 
noise  is  added  to  the  image.  Between  0  and  10  percent  noise  contamination, 
the  number  of  pixels  replaced  by  most  of  the  filters  appears  to  approach  the 
true  number  of  noisy  pixels,  but  this  is  merely  an  illusion,  for  during  the 
gradual  introduction  of  noisy  pixels,  the  respective  filters  continue  to  re¬ 
place  good  as  well  as  bad  pixels,  presumably  to  a  lesser  degree. 

Up  to  20  percent  noise  contamination,  most  of  the  filters  appear  unable 
to  keep  up  with  the  increasing  noise  contamination,  and  though  they  replace 
more  pixels  in  absolute  number,  it  is  apparent  that  a  good  deal  of  the  noisy 
pixels  are  not  being  replaced,  except  for  the  four-way  median  filter. 
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Figure  9.  Pixel  replacement  (SF11.IMG). 


From  20  to  30  percent  the  frLters  show  strong  signs  of  degradation.  The 
noise  is  beginning  to  strongly  degrade  the  bad  pixel  check  which  allows  the 
retention  of  the  noisy  pixels. 

The  four-way  median  filter  deserves  special  discussion  since  its  method 
of  application  is  unique.  In  figure  9,  the  four-way  filter  clearly  does  more 
pixel  replacement  than  any  of  the  others  at  every  percentage  of  noise  introduc¬ 
tion,  and  at  zero  percent  noise,  the  four-way  filter  replaces  more  good  pixels 
than  any  other.  As  each  pass  of  the  four-way  filter  takes  place,  fewer  and 
fewer  pixels  are  replaced.  From  the  shape  of  the  four-way  filter  envelope  in 
Figure  9,  its  general  behavior  is  similar  to  the  other  filters,  except  that  it 
performs  more  pixel  replacement,  even  at  the  higher  noise  percentages. 


The  method  of  counting  the  number  of  replacements  tor  the  four-way  filter 
is  different  than  the  method  used  for  the  other  filters.  For  the  other  fil¬ 
ters,  if  the  bad  pixel  check  indicates  that  the  pixel  needs  replacement,  the 
count  of  the  repLaced  pixels  are  incremented,  and  then  the  filtering  action  is 
taken.  For  the  four-way  filter,  if  the  bad  pixel  check  indicates  that  the 
pixel  requires  replacement  on  one  of  the  passes,  the  filtering  action  is  taken 
(via  the  five  IF  statements),  and  if  the  pixel  to  be  replaced  is  determined  to 
be  replaced  by  itself,  then  the  count  of  the  replaced  pixels  is  not  incre¬ 
mented.  The  replacement  pixel  count  is  only  incremented  when  the  pixel  to  be 
replaced  is  replaced  with  one  of  its  two  neighbors  (vertical,  horizontal,  or 
on  one  of  the  diagonals). 
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STATISTIC  FOR  PERFORMANCE  MEASUREMENT 

The  best  statistic  for  measuring  the  performance  of  the  filters  was  found 
to  be  the  standard  deviation  of  the  difference  between  the  original  image  and 
the  filtered  image  after  the  filter  had  been  applied  to  remove  the  noise. 
For  every  level  of  noise,  each  type  of  image,  whether  raw  or  scaled,  and  with 
and  without  the  bad  pixel  check,  the  correlation  coefficient  and  the  sigma  of 
the  difference  between  the  original  and  the  filtered  images  were  computed  for 
each  filter.  The  correlation  coefficient  is  not  as  good  a  comparison  statis¬ 
tic  as  the  difference  sigma.  Hence  the  difference  sigma  was  used  to  produce 
the  plots  shown  in  appendix  A. 

If  the  filtering  action  had  removed  all  of  the  induced  noise  from  all  of 
the  images,  the  mean  and  the  standard  deviation  of  the  difference  between  the 
original  and  filtered  images  would  have  been  zero,  and  the  correlation  coeffi¬ 
cient  would  have  been  1.00. 

However,  since  at  each  percentage  of  noise  replacement,  there  is  residual 
noise  after  the  filtering  action,  the  mean  of  the  difference  is  usually  approx¬ 
imately  zero,  with  a  nonzero  standard  deviation,  and  a  correlation  coetricient 
less  than  1.00.  As  increasingly  more  noise  is  placed  into  the  images,  and  as 
filtering  of  noise  decreases  as  shown  in  Figure  9,  the  mean  of  the  difference 
between  the  original  and  the  filtered  images  becomes  more  and  more  positive, 
due  to  the  greater  number  of  noisy  pixels  that  remain  in  the  filtered  image. 

FILTERING  OF  DIFFERENT  IMAGES 

With  the  exception  of  the  ocean  and  cloud  scene  (SF11.IMG  and 
SF1 1DED.IMG) ,  most  of  the  filtering  actions  upon  the  various  images  was  simi¬ 
lar.  The  four-way  median  filter  fared  better  overall  with  the  ocean  and  cloud 
scene  than  the  other  filters.  All  the  filters  had  better  performance  with  tne 
ocean  and  cloud  scene  than  the  other  images;  the  initial  disruption  of  the 
zero-noise  version  was  less,  and  from  that  point  on  with  increased  noise, 
every  filter  out-performed  itself  compared  to  its  performance  with  tne  other 
images . 

MAXIMUM  PERFORMANCE  EXPECTED  FROM  A  FILTER 

This  study  presented  an  opportunity  to  evaluate  a  given  filter  in  order 
to  determine  its  maximum  possible  filtering  capability.  Since  tne  original 
images  were  offset  from  0  by  300,  since  Gaussian  noise  was  placed  in  the  image 
in  percentages  from  0  to  30  percent,  and  since  the  value  of  noise  added  varied 
between  -20  and  +20,  a  very  simple  algorithm  was  employed  to  test  the  maximum 
potential  of  a  given  filter.  Instead  of  employing  the  usual  bad  pixel  check, 
in  its  place  a  check  was  substituted  to  determine  if  the  candidate  pixel  was 
less  than  300.  If  less  than  300,  it  was  noise  and  tne  filtering  action  was 
invoked;  if  not,  the  pixel  was  considered  good,  and  filtering  was  by-passed. 

Figure  10  compares  the  expected  and  actual  performance  of  tne  biweight 
and  neighborhood  replacement  filters.  These  two  filters  are  grouped  together 
because  they  essentially  act  the  same.  The  falloft  in  performance  as  noise  is 
increased  is  apparent  both  in  the  maximum  possible  and  actual  performance.  Of 
course,  the  actual  performance  falls  off  more  rapidly  than  the  expected.  The 
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ORIGINAL  -  FILTER  DIFFERENCE  SIGMA 

Figure  10.  Actual/expected  performance  lower  San  Francisco  Bay,  biweight  and 
neighborhood  replacement. 


closest  area  of  coincidence  is  between  2  and  4  percent  noise  contamination, 
with  the  closest  point  at  3  percent. 

Figure  1 1  compares  the  expected  and  actual  performance  of  the  median  and 
the  four-way  median  filters.  These  two  filters  are  grouped  together  because 
they  act  in  a  similar  manner.  Both  of  these  filters  perform  bad  pixel  replace¬ 
ment  with  actual  image  values,  and  not  computed  values  as  do  the  biweight  and 
the  neighborhood  replacement  filters.  The  maximum  performance  of  both  these 
filters  is  very  fine  and  very  similar.  The  actual  performance,  with  the 
exception  of  the  zero-noise  offset,  follows  the  slope  of  the  expected  perfor¬ 
mance  very  well  until  after  4-percent  noise,  when  it  degrades  in  a  manner 
similar  to  the  biweight  and  neighborhood  replacement  filters.  The  expected 
performance  of  these  filters  shown  in  Figure  11  is  mirrored  by  the  performance 
of  these  same  filters  when  the  bad  pixel  check  is  skipped  and  tne  filter  is 
applied  over  the  entire  image,  refer  to  appendix  A. 

When  the  four-way  filter  was  run  using  the  less-than-300  check  to  deter¬ 
mine  its  maximum  possible  performance,  it  was  noted  that  during  the  first 
pass,  a  greater  percentage  of  pixels  were  replaced  than  were  in  the  452x452 
image.  This  was  due  to  the  fact  that  the  area  for  the  first  pass  was  456x453 
and  thus  contained  more  noisy  pixels.  on  the  subsequent  passes,  lesser  and 
lesser  amounts  of  pixel  replacement  occurred,  indicating  that  not  all  pixel 
replacements  were  made  with  good  pixels,  but  as  the  noise  pixels  were  gradu¬ 
ally  eliminated,  the  fewer  remaining  had  a  better  chance  of  being  replaced  by 
a  good  scene  pixel. 


ORIGINAL-  FILTER  DIFFERENCE  SIGMA 


Figure  1 1.  Actual/expected  performance  lower  San  Francisco  Bay,  median  and 
four-way  median. 


IMAGE  DISTORTION  WITH  ZERO-NOISE 

The  effectiveness  of  the  various  filters  evaluated  in  this  study  depends 
also  upon  the  amount  of  distortion  caused  by  the  filter  at  zero-noise,  for  if 
any  filter  is  to  perform  well,  it  must  have  minimal  distortion  on  the  image 
when  there  is  no  noise  in  the  image.  If  there  is  a  yreat  deal  of  image  distor¬ 
tion  caused  by  the  filtering  action  when  there  is  no  noise,  then  this  tendency 
will  continue  regardless  of  the  amount  of  noise  actually  in  the  image.  An 
analysis  was  made  of  the  zero-noise  level  of  filtering  tor  all  of  the  images 
(raw  and  scaled)  and  for  the  bad  pixel  check  versus  the  no  bad  pixel  check. 
The  sigmas  of  the  difference  between  the  original  and  the  filtered  images  were 
used  as  the  criteria.  The  sigmas  for  all  of  the  filters  at  zero-noise  were 
averaged  and  were  used  to  construct  table  II.  The  average  of  the  zero-noise 
values  was  used  because  for  all  of  the  filters,  on  any  one  image,  the  sigmas 
for  each  of  the  four  filters  were  very  close  together  in  value,  as  can  be  seen 
from  the  plots  in  appendix  A.  In  most  cases  the  deviation  was  only  .01  to  .66 
in  magnitude.  The  individual  filters  were  also  ranked  by  their  aoiLity  to 
minimize  zero-noise  image  distortion. 
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Table  II.  Filter  image  distortion  with  zero-noise. 


IMAGE 

BAD  PIXEL 

CHECK 

NO  BAD 

PIXEL  CHECK 

- 

Scaled 

Sigma 

Order 

Raw 

Sigma 

Order 

Scaled 

Sigma 

Order 

Raw 

Sigma 

Order 

SF1 

3.61 

1 

2 

3=4 

1  .2 

i 

2 

3 

4 

10.8 

1  2 

3 

4 

3.59 

1  2 

4 

3 

SF2 

2.89 

4 

1 

2  3 

.97 

4 

1 

2 

3 

8.03 

1  2 

4 

3 

2.7 

1  2 

4 

3 

SF3 

7.1 

4 

1 

2  3 

2.36 

4 

1 

2 

3 

17.62 

1  2 

4 

3 

5.88 

1  2 

4 

3 

SF6 

3.75 

4 

1 

2=3 

1  .48 

4 

1 

3 

2 

9.13 

1  2 

4 

3 

3.65 

1  2 

4 

3 

SF8 

3.37 

1 

4 

3  2 

.84 

1 

4 

3 

2 

8.06 

1  2 

4 

3 

2.01 

1  2 

4 

3 

SF1  1 

1  .44 

1 

4 

2  3 

.51 

1 

2 

4 

3 

3.92 

1  3= 

4 

2 

1  .35 

1  3 

2= 

4 

Biweight=1 ,  Median=2,  Neighborhood  Replacemen t=3 ,  Four-way  Median=4 
All  sigma  values  are  averaged  for  all  four  filters. 

One  outstanding  characteristic  of  all  of  these  filters  is  that  they  all 
have  the  same  propensity  of  image  distortion  at  zero-noise.  Even  though  the 
biweight  filter  was  better  than  the  rest,  its  performance  in  terms  of  the 
others  was  not  that  much  better;  the  ranking  is  only  an  artifact  of  comparing 
the  absolute  values,  and  being  forced  to  select  one  as  greater  than  or  less 
than  another.  In  view  of  this  comment,  the  performance  of  the  filters  at  zero- 
noise  is  (1)  biweight,  (2)  median,  (3)  four-way  median,  and  (4)  neighborhood 
replacement,  from  best  to  worst. 

At  zero-noise,  the  effect  of  not  performing  the  bad  pixel  check  is  to 
introduce  more  filter  distortion  into  the  image,  however  this  tendency  is  not 
so  pronounced  when  comparing  the  raw  images.  At  zero-noise  approximately  3 
percent  of  the  pixels  are  distorted  from  the  original  and  if  as  mentioned 
earlier,  some  of  this  distortion  is  actually  removing  noise  that  was  present 

in  the  original  image,  then  this  concern  of  distortion  at  zero-noise  would 

become  lessened.  However,  when  the  filter  is  applied  to  the  entire  image 
(without  the  bad  pixel  check),  the  median  and  four-way  median  filters  do  tend 

to  follow  their  maximum  possible  filtering  curve  better  than  when  the  bad 

pixel  check  is  performed. 

IMPROVED  BAD  PIXEL  CHECK 

The  bad  pixel  check  was  improved  by  increasing  the  neighborhood  to  11x11 
and  using  various  number  of  sigmas  for  increasing  amounts  of  induced  noise  as 
shown  in  table  III. 
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Table  III.  Improved  bad  pixel  check  (11x11  neighborhood) 


Noise  (%) 

#  of  Sigmas 

Noise  Pixels 

Noise  Replaced 

Good  Disturbed 

0 

4.1 

0 

0 

14 

1 

3.5 

2076 

2076 

24 

2 

3.5 

4074 

4074 

14 

3 

3.5 

6021 

6021 

6 

4 

3.1 

8035 

8035 

1  9 

5 

3.0 

10,069 

10,068 

1 1 

6 

2.9 

12,084 

12,079 

5 

7 

2.8 

14,056 

14,047 

5 

8 

2.7 

16,068 

16,061 

3 

9 

2.6 

18,075 

18,067 

3 

10 

2.5 

20,110 

20,105 

5 

20 

1  .5 

40,259 

40,259 

0 

30 

1  .3 

60,716 

60,712 

0 

40 

1  .0 

81 ,424 

81 ,424 

16 

Tests  were  run  on  the  raw  DAEDALUS  images  downtown  San  Jose  (SF3DED)  and 
ocean  and  clouds  off  Santa  Cruz  (SF11DED).  These  are  respectively  the  most 
and  least  cluttered  of  the  images.  The  results  are  shown  in  figure  12  for 
downtown  San  Jose,  and  in  figure  13  for  ocean  and  clouds  off  ->anta  Cruz.  The 
improvement  in  both  images  over  the  performance  using  the  bad  pixel  check  with 
a  5x5  neighborhood  and  constant  number  of  sigmas  =  2  is  substantial.  The 
unevenness  of  the  plots  between  5-  and  10-percent  noise  could  be  eliminated  by 
fine-tuning  the  selection  of  number  of  sigmas,  because,  as  table  HI  shows, 
not  all  of  the  noise  was  removed.  Residual  noise  contributes  more  to  in¬ 
creasing  the  difference  sigmas  (difference  between  the  original  and  filtered 
images)  than  does  the  small  amount  of  good  pixels  disturbed. 


Figure  12.  Improved  filter  performance  (SF3DED.IMG). 


SUMMARY  OF  FILTER  PERFORMANCE 


The  major  consideration  in  arriving  at  a  final  determination  of  filter 
performance  is  not  to  look  at  the  statistics  alone,  but  also  to  look  at  the 
final  filtered  image.  All  of  the  filters  have  a  tendency  to  blur  or  distort 
the  original  in  some  manner.  The  manner  of  inducing  noise  in  this  study  was 
not  to  add  or  subtract  some  constant  (or  even  variable)  to  the  pixels  in  the 
image.  The  actual  pixel  values  were  replaced  with  values  that  have  no  rele¬ 
vance  to  the  original  pixel  value.  The  original  pixel  value  is  gone;  it  can 
never  be  completely  recovered.  The  filtering  techniques  were  evaluated  to  see 
how  well  the  substituted  value  came  close  to  the  original  pixel  value. 

Another  consideration  is  the  amount  of  noise  that  a  filter  is  expected  to 
replace.  Even  with  levels  of  noise  as  high  as  50  percent,  it  is  surprising 
that  a  human  being  can  still  recognize  the  major  content  of  an  image.  In  this 
respect,  all  of  the  filtering  techniques  appear  to  be  performing  well.  On  the 
other  hand,  if  restoration  of  an  image  is  required  to  make  intelligible  some 
small  2x2  or  4x4  pixel  area  for  identification  purposes  when  50  percent  of  the 
pixels  are  replaced  by  noise,  then  these  types  of  filters  would  be  inadequate, 
if  indeed  any  would  be  adequate. 

Several  criteria  have  been  selected  to  evaluate  the  filters  in  view  of 
the  above: 

1.  The  potential  of  the  filter  at  different  levels  of  noise 

2.  The  effectiveness  of  the  filter  at  different  levels  of  noise 

FILTER  POTENTIAL 

Figures  10  and  11  show  that  the  maximum  potential  for  the  tour  filters  is 
limited  assuming  that  a  bad  pixel  check  can  be  designed  which  will  precisely 
identify  only  the  bad  pixels,  and  if  these  are  compared,  the  ranking  of  the 
filters  would  be  (1)  median,  (2)  four-way  median,  (3)  biweight,  and  (4) 
neighborhood  replacement.  This  is  the  best  that  the  filters  could  perform. 
The  expected  performance  is  thus  the  best  it  could  possibly  be  throughout  the 
entire  range  of  noise  replacement,  from  0  percent  to  30  percent. 

Viewing  these  plots  shows  that  the  consistency  of  maintaining  a  good 
replacement  of  noise  with  values  close  to  perhaps  the  original  pixel  values  is 
continuous  throughout  the  range  for  the  median  and  four-way  median  filters.  On 
the  other  hand,  the  biweight  and  the  neighborhood  replacement  filters  fall  off 
in  potential  performance  as  the  amount  of  noise  contamination  increases.  The 
median  and  the  four-way  median  filters  replace  noise  with  actual  original 
pixel  values,  whereas  the  biweight  and  the  neighborhood  replacement  filters 
manufacture  new  pixel  values  that  may  never  have  existed  in  the  original  image. 

FILTER  EFFECTIVENESS 

Until  the  bad  pixel  check  was  improved  to  reduce  or  eliminate  zero-noise 
image  distortion,  the  performance  of  all  four  filters  throughout  the  entire 
range  of  0  to  10  percent  noise  fell  short  of  the  expected  performance,  as 
shown  by  Figures  10  and  11.  The  four-way  median  and  median  filters  performed 


well  in  the  0-  to  4-percent  noise  range,  so  they  could  be  classed  as  the 
better  filter. 


Without  the  bad  pixel  check,  clearly  the  median  and  four-way  filters  are 
more  effective  than  the  biweight  and  neighborhood  replacement  filters  as  shown 
in  figure  8. 

Figures  12  and  13  show  that  the  median  and  four-way  median  filters  out¬ 
perform  the  biweight  and  the  neighborhood  replacement  filters,  and  that  the 
four-way  median  filter  outperforms  the  median  filter  at  noise  levels  higher 
than  10  percent.  These  plots  were  made  using  the  improved  bad  pixel  check 
discussed  earlier. 

SURVEY  OF  FILTERED  IMAGES 

A  program  was  written  to  display  the  original  image,  the  original  image 
contaminated  by  noise,  the  noise  mask,  the  filtered  image,  and  the  mask  of  the 
residual  noise  and  good  pixels  disturbed  by  filtering.  Due  to  the  fact  that 
in  many  cases  residual  noise  remained  in  the  filtered  image,  the  image  would 
have  to  be  rescaled  in  order  to  display  it.  The  rescaling  with  noise  values 
from  -20  to  +20  distorted  the  filtered  image  compared  to  the  original.  It  was 
decided  to  produce  a  synthetic  representation  of  the  filtering  action.  The 
original  image  was  input  from  disk  and  displayed  as  is.  Then  another  working 
image  was  built  of  the  original  image  plus  the  300  offset  and  plus  noise. 
This  was  the  image  upon  which  the  filter  worked.  The  original  image  had  the 
same  noise  pixels  as  the  working  image,  but  instead  of  values  from  -20  to  +20 
the  values  were  what  the  operator  selected,  usually  zero,  which  displayed  as 
black.  The  original  image  with  black  noise  was  displayed,  then  a  noise  mask 
was  displayed.  The  noise  mask  was  built  at  the  same  time  the  black  noise  was 
added  to  the  original  image.  The  noise  mask  was  black  with  white  noise.  Then 
the  filtering  action  took  place.  As  each  pixel  requiring  replacement  (as 
determined  by  the  bad  pixel  check)  was  corrected  in  the  working  image,  a  check 
was  made  in  the  original  image  for  presence  of  zero  pixel  values  (indicating 
noise).  If  the  pixel  to  be  replaced  was  zero,  then  the  count  of  noise 
replacement  was  incremented  and  the  noise  pixel  in  the  mask  was  zeroed;  if  not 
zero  the  count  of  pixels  disturbed  was  incremented  and  the  pixel  in  the  noise 
mask  was  made  purple. 

In  either  event  the  pixel  in  the  original  image  was  replaced  with  the 
replacement  value  in  the  working  image  minus  300.  This  was  not  the  correct 
value,  but  served  to  indicate  that  a  correction  had  been  made,  and  also  re¬ 
moved  the  zero  (noise)  values  from  the  original.  After  filtering,  the  origi¬ 
nal  with  synthetic  noise  removal  was  displayed,  and  finally  the  residual  noise 
mask  was  displayed.  Pictures  were  taken  of  the  sequences  for  different  images 
and  varying  amounts  of  noise.  The  filter  used  for  the  filtering  action  was 
the  median  filter. 

Figures  14,  15,  and  16  show  respectively  1-,  5-,  and  50-percent  noise 
masks.  It  is  not  likely  that  any  images  containing  50-percent  noise  would  be 
processed,  but  figure  16  gives  an  idea  of  that  amount  of  noise. 


Figure  17  shows  1 -percent  noise  applied  to  the  ocean  and  cloud  image  off 
Santa  Cruz.  The  picture  was  dark  so  white  noise  was  synthetically  added. 
Below  and  to  the  right  of  the  image  center  is  a  ship.  Noise  is  not  added  to 
the  edge  of  the  image. 

Figure  18  shows  the  results  of  filtering  using  a  3x3  neighborhood  and 

number  of  sigmas  =  2  for  the  bad  pixel  check.  Original  noise  remains  in  the 
five-pixel  border,  which  the  filter  does  not  cover,  and  some  residual  noise  re¬ 
mains  as  seen  by  the  white  noise  still  in  the  image. 

Figure  19  shows  the  residual  noise  mask.  The  large  dots  are  the  residual 
noise;  the  small  dots  are  good  pixels  which  have  been  disturbed.  In  black  and 
white  photos  the  difference  between  the  residual  noise  and  the  good  pixels 
that  have  been  disturbed  does  not  stand  out  too  well. 

Figure  20  shows  5-percent  noise  replacement  in  the  lower  San  Francisco 

Bay  image.  The  noise  is  in  black. 

Figure  21  shows  the  results  of  median  filtering.  The  bad  pixel  check 

employed  a  9x9  neighborhood  and  number  of  sigmas  =  2.5.  Of  the  10,069  noisy 
pixels  in  the  original,  the  filter  removed  10,040,  and  disturbed  only  24  good 
pixels.  The  residual  noise  mask  for  this  filtering  is  shown  in  figure  22. 

Figure  23  shows  northern  San  Jose  with  10  percent  noisy  pixels  in  black. 
Using  a  9x9  neighborhood  and  number  of  sigmas  =  2.0,  19,996  noisy  pixels  were 

replaced  out  of  20,110,  and  no  good  pixels  were  disturbed.  The  results  of 
filtering  are  shown  in  Figure  24  and  the  residual  noise  mask  is  shown  in 
figure  25. 

Figure  26  shows  Santa  Cruz  and  tne  ocean  with  30  percent  black  noise 

added.  Using  9x9  neighborhood  and  number  of  sigmas  =  1  .0,  60,709  of  the 

60,719  noise  pixels  were  replaced  by  the  median  filter,  and  only  two  good 
pixels  were  disturbed.  The  results  of  the  filtering  is  shown  in  figure  27, 
and  the  residual  noise  mask  is  shown  in  Figure  28. 


ff  Santa  Cruz  median  filtered. 


Figure  25.  Residual  noise  on  northern  San  Jose 
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RECOMMENDATION  FOR  FUTURE  STUDY 


Another  improvement  to  the  bad  pixel  check  would  be  to  build  a  bad  pixel 
mask  as  shown  in  figure  29.  This  would  be  done  by  applying  the  bad  pixel 
check  to  the  entire  image  three  times.  On  the  first  pass  the  bad  pixel  check 
would  sum  only  the  good  pixels  in  its  neighborhood  for  computation  of  the  mean 
and  sigma.  Since  the  first  time  through,  the  pixels  above  and  to  the  left  of 
the  candidate  would  have  already  been  marked  as  a  bad  or  good  pixel,  the  bad 
pixels  can  be  left  out  of  the  sum. 
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BAD  PIXEL  MASK 


FIRST  PASS: 

Compute  mean  and  sigma  of 
submatrix  using  only  good  pixels 
mark  bad  pixels  as  1 . 


SECOND  PASS: 

Skip  good  pixels,  recheck  bad 
pixels  using  only  good  pixels. 
Make  bad  pixels  good  if  within 
selected  number  of  sigmas 
from  mean. 


THIRD  PASS: 

Use  bad  pixel  mask  to  filter 
and  replace  only  bad  pixels 
(as  per  maskl. 


Figure  29.  Bad  pixel  mask. 


On  the  second  pass,  all  bad  pixels  in  the  mask  can  be  checked  to  see  if 
they  are  really  bad  pixels;  if  not,  they  can  be  returned  to  the  status  of  good 
pixels.  The  mean  and  sigma  computations  would  be  made  using  only  the  good 
pixels  in  the  selected  neighborhood.  On  the  third  pass,  actual  filtering  can 
be  carried  out  using  a  3x3  neighborhood,  but  only  the  good  pixels  in  that 
neighborhood  would  be  processed  to  replace  the  marked  bad  pixels.  On  the 
third  pass,  only  the  bad  pixel  mask  would  be  screened  to  determine  which 
pixels  require  replacement.  Time  would  be  saved  on  the  third  pass,  since  the 
noisy  pixels  have  already  been  tagged.  The  filters  would  use  only  the  good 
pixels  in  the  3x3  neighborhood  around  the  candidate  pixel.  Variations  in 
filtering  would  have  to  be  made  for  each  filtering  operation.  The  biweight 
filter  would  have  to  be  revised  to  operate  on  less  than  nine  pixels.  The  sort 
routine  SHELLSORT  is  designed  to  handle  any  number  of  values. 

The  median  filter  would  have  to  be  revised  to  handle  even  numbers  of 
pixels  (2,  4,  6,  or  8)  by  averaging  the  first  and  second  (for  n=2 ) ,  the  second 
and  third  (for  n=4),  the  third  and  fourth  (for  n=6),  and  the  fourth  and  fifth 
(for  n=8);  and  using  this  average  for  replacement. 


The  neighborhood  replacement  filter  would  sum  only  the  good  pixels  in  its 
neighborhood  and  divide  by  the  n  of  the  good  pixels. 

The  four-way  median  filter,  on  each  of  its  passes,  would  not  replace 
unless  at  least  one  good  pixel  was  found  in  one  of  the  two  wings  of  the  ver¬ 
tical,  horizontal,  or  diagonal  sets. 


In  all  the  filters,  when  a  replacement  is  made,  the  bad  pixel  mask  is 
updated  by  making  the  bad  pixel  good.  This  means  that  the  replaced  pixel  can 
be  used  subsequently  in  the  computations  for  replacing  other  bad  pixels  in  its 
3x3  neighborhood. 
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program  onfl 

this  program  roads  daedalus  images  from  disk  than  introduces 
normal  noise  then  applies  biueight  filter  to  remove  the  noise 
then  computes  correlation  between  original  and  filtered  images 
include  ' teamdef.  f or/nol i st  ' 

INTEGER*2  cc(262144)  ‘corrupted  image 
integer»4  seed,  r  012) 

real* S  dsum.  dsq.  corr.  so.  sc.  poc.  ssqo.  ssqc 
real  med<9>  ! array  for  sorting 
bgte  dd(262144)  ! original  image 
character*20  f  .'file  name 
seed =77737 
do  1  i=l. 512 

r  (  i  )  =  (  i-1  )*512  .'build  row  table 

i=GT_INITIAL< 1. 2)  ! init  term  function  **************************** 

call  GT_PRQMPT< '0NF1>  ')  ! terminal  prompt 

i=gt_word( 'Disk  File  Name < SFn(DED) .  IMG >  =  ' .  f < 1 :  20) > 

open<unit=7, file=f. form* 'unformatted'. status='old'>  err =88) 

read (7, err=B8). dd  ! GET  IMAGE 

type  *.  'IMAGE=  '.  f 

c lose ( uni t=7) 

do  7  ir=26. 487  (converting  image  from  bgte  to  i*2 
do  7  ic=26. 487 

ip-dd (r ( ir )+ic >  forig  bgte 

if(ip  .  It.  0) ip=ip+256  !  if  128-255  make  positive 
ip=ip+300  (shift  up  to  get  awag  from  0  mean  gausian  noise 
cc (r ( ir )+ic )=ip  (put  in  contaminated  image. tho  not  get  contam 
i=gt_real('#  Sigmas  for  bad  pixel  cr iterion* '. sig ) 
i=gt_integ ( 'Neighb orhood  Size  <ie  3x3=1 . 5« 5=2 >=', is ) 
i=gt_integ < '0=Bad  Pix  Check.  l=Skip  ='.ik> 
fn=< is+i s+1 >*< is+i s+1 )  'compute  n  for 
g =204304. 0  (#  of  pixels  in  452*452  image 
ibp=0  ! #  pixels  contaminated 

if <gt_ask ( 'Want  to  introduce  noi se '. none >  eq. no > goto  67 
i“9*_r*al < 'Enter  Percent  contaminated '. p > 
i“9t_real < 'Enter  Gaussian  Mean'.gm) 
i=gt_real ( 'Gaussian  std.  dev.  ',sd> 

do  66.  ir=26. 487  (  this  loop  enter  1-p  7.  gaussian  noise  mean  O.  sigma=sd 

do  66.  ic=26. 487  (corrupt  whole  picture 

if(ran(seed). gt. 1 . 0-p ) then 

cc<r< ir)+ic)=nint<gauss(gm. sd) )  ! noise 
if<ir  .  ge.  31  and.  ir  Ie.  482>then 

if<ic  .  ge  31  and  ic  le.  482) ibp=ibp+l  !#  pixels  corrupted 
end  i  f 
end  i  f 
continue 
dsum=0  O 
dsq=0.  O 
do  8  ir=31. 482 
do  8  ic=31« 482 
ip=dd<r<ir)+ic) 
iflip  It  0)ip=ip+256 

pd=f 1  oat < ip+300 )-f 1 oa t ( cc < r < i r > +ic > >  (dif  orig-contam 
dsum=dsum+pd  (Sum  of  diffs 
dsq=dsq+pd*pd  'sum  of  diff  squared 

dm*an=dsum/g  (mean  of  difference  of  orig  and  contaminated  images 

d sd=sqr t ( < d sq/g ) -dmean*dmean >  !  stan  dev  of  difference 

type  *. 'Diff(CONTAM  -ORIG  )  Mean*'. dmean 

type  *.  ' Di f f ( CONTAM  -ORIG  )  Sigma='.dsd 

icp=0  (W  pixels  corrected 

igp=0  good  pixs  disturbed 


nbp90  1  #  bad  pixs  corrected 

do  25, ir=31.482  ! detect  bad  pixels  and  correct 
do  25  .  ic=31.  482 

iflik  ne.  O)goto  23  (skip  had  pixel  check,  filter  all 
i=float(cc(r< ir )+ic ) ) 

dsq«0  0 
dsum=0  0 

do  21  iro=-is. is  1 fn  window  compute  mean  and  sigma 
do  21  ico“~is, is 

pc =f loat ( cc < r ( i r  +  iro) +i c  +  i c o ) )  !  get  pixel 
dsum“dsum+pc  1  sum 
dsq«dsq+pc*p c  f sum  squared 
cmean-dsum/f n  !  compute  local  mean 
csd  =  sqr  t  <  (ds  q/fn  )-'cmean*cmean  1  Moca)  stan  dev 

if  <  <  x  ge  cmean-s i g*c  sd 1 .  and.  <  x  1  e  cmean+sig*csd))goto  25  !  quod 
icp=icp+l  ■ count  corrected  pixels  subtract  from  contaminated 
ict=0 

i f < c c < r < ir >♦  ic  1  It.  200)then 
nbp=nbp+l  (had  pixs  corrected 
else 

igP-i9P+l  ‘good  pixs  disturbed 
end  i  f 

do  22  iro=-l.l  (extract  3x3  window 
do  22  ic o=~ 1 .  1 

p c “float < cc < r < i r+i ro )+ l c +i c o > )  'get  pixel 
l  c  t  =  ic  t+ 1 

med(ict>=pc  ( med  index  varies  from  1-9 
call  sh e  1 1  so r t ( med .  9 )  1  veq  for  biugt  .V  ined 
call  b iweigh t (med.  x > 
cc(r(ir)+ic)“nint(x> 
continue 

type  *.  '#  noise  pixels9',  ibp.  '  #  corrected  pixels-',  icp 
type  *,  '#  noise  pixels  fixed9', nbp,  '  #  good  disturbed9',  lqp 
so=0. 0  (Sum  of  orig 
sc=0  0  (Sum  of  filtered 
poc“0.  O  1  cross  prod 
ssqo=0. 0  (sum  orig  squared 
ssqc=0. O  (sum  filt  squared 
dsum-O.  O  'sum  of  diff 
dsq=0  O  'sum  dif  squared 
do  18  ir=31 . 482 
do  18  ic-31.482 

i p=dd <r < ir )+ic )  ‘original  pixel 
if (ip  It  O  > i p  =  i p+25A 
po“float< ip +300) 
pc  =  float(cc(r< i r ) +i c )  ) 
so=so+po  '  sum  of  orig 
sc“sc+pc  'sum  of  filtered 
poc “poc +p o*pc  (sum  cross  product 
ssqo“ssqo+po*po  (  sum  orig  squared 
ssqc»ssqc+pc*p c  'sum  filtered  squared 
dsum“dsum+(po-pc )  (sum  of  diffs 
dsq=dsq+< po-pc )*( p o-pc )  (sum  of  dif  squared 

dmean“dsum/g  (mean  of  difference  of  original  and  filtered  images 

dsd=sqr t ( ( dsq/g )-dmean*dmean )  ! stan  dev  of  difference 

corr  =  (g*poc-so*sc)/<sqrt<  <  g*ssqo-so*so )*<  g*ssqc-sc*sc ) 1 ) 

type  *.  'CORRELATION) or i g inal .  f i 1  ter ed )  =  c orr 

type  *. 'Di f f < F ILTERED-OR IG >  Mean= ' . dmean 

type  *.  'Di f f (FILTERED-OR IG 1  iii gma=  d  s d.  '  * 

type  *.  'Bi weight  FILTER  NOISE9  p 

type  *.  'NOISE  Clean  =  '.  gm.  '  NOISE  Std  de v=  ' .  sd ,  '  Ne  i  gh  -  ' .  is 

call  exit 
end 

include  'mathfunc  for/nolist' 
include  'shellsort  for/nolist' 
include  'biweight  for/nolist' 


p»v»«n 


sew 


wrap 


program  onif' 

c  this  program  reads  daedalus  image*  from  disk  then  introduces 

c  normal  noise  then  applies  the  median  filter  to  remove  the  muse 

c  then  computes  correlation  between  original  and  filtered  images 

c  option  added  to  skip  bad  pixel  check 

c  variable  neighborhood  for  bad  pixel  check  m  3x3.  5x5.  etc 

include  'teamdef  for/nolist' 

INTEGER*2  cc<262144)  'this  is  corrupted  image 
integer*4  seed.r<512> 

real*8  dsum.  dsq,  covri  so.  sc*  poc,  ssqo,  ssqc 
real  med(9)  1  array  for  sorting 
byte  dd(262144)  'original  image 
character*20  f  ‘file  name 

seed  =77  737 
do  1  i=l,512 

1  r ( i )  =  < i- 1 ) *5 12  !  hoild  row  table 

c  *##■####*#####*###«  ***#*##**•##*###***•*#■*#***#*•*#*#*##*#**#■##*##** 

i=GT_INITIAI  (1,2)  1  init  term  function 

call  GT_PRuMPT (  '0NF2>  ')  1  terminal  prompt 

c  ******  ************  ********  ******************************* ******* 

88  i=g t_word < 'Disk  File  Name ( GFn ( DED )  IMG )  =  ' ,  f < 1  20 )) 

open  (  uni  t =7*  f  i  le-=f ,  torm-  'unformatted  status-  old  .  er  r=80 ) 
read ( 7, err=88) , dd  ‘GET  IMAGE 
type  *.  '  IMAGE*  ',  f 
close  <unit*7) 

do  7  ir-26, 487  'converting  image  from  byte  to  i*2 
do  7  ic=2fc, 487 

i  p  ~  d  d  (  r  (  i  r  )  + 1  c  )  1  ong  byte 

if(ip  It  0)ip  ip+256  1  if  128-255  make  positive 
ip=ip+3G0  'shift  up  to  get  away  from  O  mean  gausian  noise 
7  cc  (rur)  *u  >  =  ip  'put  in  contaminated  image,  tho  not  yet  contain 

i=gt_real<#  Sigmas  for  bad  pixel  criterion”',  sig) 
i=gt_integ< 'Neighborhood  Size  (le  3x3=1, 5x5-2)=", is) 
i  =  gt_integ(  ‘0--Bad  Pix  Check;  l=Skip  =',  ik> 
fn=( is+15+1 )#< is+i s*l )  ! compute  n  for  bad  pix  check 
g=204304  O  '  #  of  pixels  in  452*452  image 
lbp-0  '  #  pixels  corrupted 

i  f  (  g  t  _as  k  (  Wan  t  to  introduce  no  l  se  7 ,  none  )  eq  no)guto  67 
i=g t_real (  'Enter  Percent  t on tamina t ed ' ,  p ) 
l  s9  t..rea  1  (  'Gauss  i  an  means',  gm)  'allow  any  gaussian  mean 
i  =  gt_.real  (  'Gaussian  std.dev.  ',sd) 


do  66,  ir=26»  487  !  this  loop  enter  1-p  V.  gaussian  noise  mean  0,  siq,na=sd 

do  66.  ic  -  26,  487 

l f  < ran (seed),  gt  1  Op)then 

cc  <r(  i r ) *i c  )=nint<gauss(gm,  sd)  )  '  noise 
if  (  ir  ge  31  and  it-  le  482)then 

if ( ic  ge  31  and  ic  le  482 > i b p  - 1 b p  + 1  '4  pixels  corrupted 
end  if 
end  i  f 
continue 


t************** 


dsum=0  0 
d sq-0  0 

do  8  i r =31 , 482 
do  8  1 c  =3 1 , 402 
ip  =  dd(r(ir)Mc  ) 
lfCip  It  0)  l  p  -i  p-»256 


pd=f loat < ip+300)-f lost (cc < r ( ir )  +  i c ) )  Idiff  orig-corrupted 
dsumadsum+pd 
S  dsq=dsq+pd*pd 

dmean=dsum/g  ! mean  of  difference  of  orig.  and  contaminated  images 

dsd=*sqr t < ( dsq/g )-dmean*dmean )  ! stan.  dev.  of  difference 

type  *, 'Dif f (CONTAM. -ORIG.  >  Mean='.dmean 

type  *, 'Dif f (CONTAM. -ORIG  )  Sigma=',dsd 

icp-0  !  #  pixels  corrected 

nbp“0  t 

igp”0 

c ******************************************************************* 
do  25  ir«31, 482  !  detect  bad  pixels  and  correct 
do  25  i c -31 >  4B2 

if<ik  . ne.  O)goto  23  !  skip  bad  pixel  check/  filter  all 
x“float<cc(T<irl-*-ic)  l 
dsq^O  0 
dsom=0.  0 

do  21  iro«-is. is  ! fn  window  compute  mean  and  sigma 
do  21  ico=-is. is 

pc=f loat ( cc ( r ( ir+ir o >+ic+ic o > )  !  get  pixel 
dsum*dsum+pc  (  sum 

21  dsq=dsq+pc*p c  ‘sum  squared 
cmean-dsum/f n  ! compute  local  mean 

csd=sqr t (dsq/fn-cmean*cmean >  !  local  stan.  dev 
i  f  (  (  x  ge.  c mean- sig*csd>.  and.  (x.  le. cmean+sig»csd ) ) goto  25  !g 
23  icp=icp+l  .'count  corrected  pixels 

if ( x  It.  100  0) then 

nbpanbp+1  ! found  bad  pixel 
else 

igpaigp+1  ‘ correcting  good  pixel 
end  i  f 
ict*0 

do  22  iro»-l.l  ! 3x3  window  filter 
do  22  ico=-l.  1 

pc=f loat < cc <r < ir+iro >+ic+ico ) )  !  get  pixel 
ict=ict+l 

22  med(ict)=pc  ! med  index  varies  from  1-9 
call  shel lsort(med/ 9)  ! req  for  med 

cc  ( r  <  ir  ) +  ic  >  =*nint  (med  (  5)  >  1  median 

i  f  (  f I  oat ( cc ( r ( ir  )  +  i  c >  ) .  e  q  x)igp=igp-l  ’  same  pixel 
25  continue 

C ************************************* ****************************** 
type  *,  '*  corrupted  pixels*",  ibp/  '  bcorrected  pixels*".  icp 
type  */  '#  corrected  p i xe 1 s= ' . nbp.  '  tdisturbed  pixels-'/  igp 
>0*0  O  ’Sum  of  orig 
sc«0.  O  ! Sum  of  filtered 
poc=*0.  0  !  cross  prod 
ssqoaO.  O  !  sum  orig  squared 
ssqc**0.  0  (sum  filt  squared 
dsum=0. O  !  sum  of  diff 
diq*0  0  ‘sum  dif  squared 
do  18  ir»31,482 
do  18  ic=31 . 482 

ip>dd ( r ( i r ) +i c )  ioriginal  pixel 
if(ip  It  0) i p=ip+256 
po=float( ip *300) 
pc“float(cc (r< ir)+ic ) ) 
soaso+po  ( sum  of  orig 
sc=sc+pc  'sum  of  filtered 
pocapocrpoepc  'sum  cross  product 
ssqo*ssqo+po*p o  'sum  orig  squared 
ssqc=ssqc +pc*p c  'sum  filtered  squared 
d sum*d sum+ < po-p c )  'sum  of  diffs 
18  dsq»dsq-Mpo-pc )*(po-pc )  'sum  of  dif  squared 

dmeanadsum/g  'mean  of  difference  of  original  and  filtered  images 
dsd«sqrt << dsq/g >-dmeen*dmean )  'stan  dev  of  difference 
corra(g*poc-so*sc )/<sqrt<  (g*ssqo-so*so)*(g*ssqc-sc»sc  )  )  ) 
type  •/  'CORRELATION or l g ina 1 .  f 1 1  ter ed )  =  ' , c orr 
type  *. 'Di f f (FILTERED-ORIGl  Meen='.dmean 

type  *.  'Diff (FILTERED-ORIO  Sigma='.dsd,  » 

type  e, 'Median  FILTER  NOISE*". p 

type  *. 'Gaussian  NOISE  Mean='.gm. '  NOISE  Std  Dev='.sd 

call  exit 

end 

include  'mathfunc  for/nolist' 
include  'shellsort  for/nolist' 

C-5 


n  ni 


program  onf3 

this  program  reads  daedalus  images  from  disk  then  introduces 
normal  noise  then  applies  the  neighborhood  replacement  filter  to 
remove  the  noise  then  computes  correlation  between  original  and 
filtered  images 
include  ' teamdef .  f or/nol i st ' 

INTEGER*2  cc<262144)  'this  is  corrupted  image 
real*8  dsum.  dsq,  corr.  so.  sc.  poc.  ssqo,  ssqc 
integer*4  seed.r(512) 
byte  dd<262144>  'original  image 
character*20  f  'file  name 


seed=77737 
do  1  i  =  l, 512 

1  r ( i )=( i-1 )*512  'build  row  table 

C ************************* ************************************ 

i=GT_INITIAL < 1 .  2)  !init  term  function 

call  GT_PROMPT< '0NF3>  ')  (terminal  prompt 

C ************************************** *********************** 

88  i=gt_word<  'Disk  File  Name (SFn < DED > . IMG >= '. f < 1  20 > > 

open(unit=7, fil e=f . f orm= 'unformatted', status='old'. prr =88 ) 
read<7, err=88>. dd  'GET  IMAGE 
type  *,  '  IMAGE=  '.  f 
c lose ( unit-7) 

do  7  ir=26. 487  'converting  image  from  byte  to  i*2 
do  7  ic=26. 487 

ip=dd <r < ir )  +  i c  )  !orig  byte 

i f  < i p  .  It.  0)ip=ip+256  '  if  128-255  make  positive 
ip=ip+300  ‘shift  up  to  get  away  from  O  mean  gausian  noise 
cc ( r ( ir )  +  ic )=ip  'put  in  contaminated  image. tho  not  yet  contam 
************************************************************ 
i=gt_r*al('#  Sigmas  for  bad  pixel  c r i ter i on=  ' , s i g > 


i=gt_ 

integ  <  'Neighb 

orhood 

Si  le 

< ie  3x3=1, 5x 

5=2)=' 

/is) 

i*gt_ 

integ ( '0=  Bad 

Pix  Check, 

l=Sk ip ' .  i k  ) 

fn=(  i 

s+is+1 )*< is+i 

s+1)  ' 

n  for 

neighborhood 

g=204304  0  ' «  of  p 

i  xe  1  s 

in  452*452  image 

ibp=0 

' #  pixels  co 

rrup  ted 

i  f  ( g  t 

_ask ( 'Want  to 

intro 

duce  noise', none) 

eq  no ) 

goto  67 

i  =g  t_ 

real < 'Enter  P 

ercent 

contaminated  ', p  ) 

i=gt_ 

real ( 'Gaussia 

n  mean 

=  ' *  gm) 

1  allow  any 

gaussi 

an  mean 

66 

C*« 

67 


i=gt_real < 'Gaussian  std  dev  '.sd) 

»****»*********contaminate  with  noise*********************** 
do  66.  ir=26. 487  '  this  loop  enter  1-p  V.  gaussian  noise  mean  0. sigma=sd 

do  66.  ic=26,  487 

if (ran(seed)  gt  1  0-p)then 

cc(r(ir)+ic)=nint(gauss<gm,sd))  ' noise 
if(ir  go.  31  and.  ir  le.  482)then 

if<ic  ge  31  and  ic  le  482 )ibp=ibp+l  '#  pixels  corrupted 
end  i  f 
end  if 
continue 

******* *********d if  between  on g /c ontami na t ed ******************** 
dsum-0.  0 
dsq=0  O 

do  8  l r=31 .  482 
do  8  l c=31 . 482 
ip=dd<r(ir)*ic ) 
if(ip  It  0)ip=ip+256 

pd  =  f 1  oat < i p+300 )-f 1  oat < c c ( r < i r ) *ic ) )  'diff  or  l  g-c or r up t ed 


d  sum=d  sum+p  d 


C-6 
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dsq=dsq+pd*pd 

dmean=dsum/g  ! nun  of  difference  of  orig.  and  contaminated  images 
dsd*sqrt ( < dsq/g )-dmean*dmean )  ! stan.  dev.  of  difference 
type  *.  'Di  f  f  (CONTAM.  -ORIG.  )Mean=  ' »  dmean 
type  *. 'Diff <CONTAM  -ORIG  )Sigma=',dsd 
icp»0  ! H  pixels  corrected 
igpsO  !#  good  pixels  disturbed 
nbpsO  ! ttnoise  pixels  corrected 
c************»**FILTER  IMAGE  ************************************** 
do  25  ir=31,4B2  ! detect  bad  pixels  and  correct 
do  25  ic=3t< 482 

x=float<cc<r(ir>*ic>> 

i f  < i k  .  ne  Olgoto  23  Iskip  bad  pixel  check 
dsq=0  0 
dsum=0  0 

do  21  iro--isiis  ! 11x11  window  compute  mean  and  sigma  for  window 
do  21  ico=-is> is 

pc=f loat < cc < r < ir+ir o > *i c+i co > >  ! get  pixel 
dsum=dsum+pc  ! sum 

21  dsq=dsq+pc*pc  ! sum  squared 
cmean=dsum/f n  !  compute  local  mean 

c sd=sqr t ( d sq/f n-cmean*cmean )  ! local  stan. dev 

i  f  <  <  x.  ge. cmean-s ig*c  sd ) .  and  <  x  lecmean+sig*csd))goto  25  !  good 
23  icp=icp+l  !  count  corrected  pixels 

if <x  It.  100  O) then 

nbp=nbp+l  ! noise  corrected 
else 

igp=igp*l  'good  disturbed 
end  i  f 
dsum=0  0 

do  22  iro=-l.l  ! 3x3  window  filter 
do  22  ico=-l.  1 

22  dsum=dsum+pc  summing  pixels  1-9 

cc < r ( i r ) +ic ) =nint < < dsum-x > /8. 0)  ‘  neighborhood  replacement 
i f ( f 1  oat (cc ( r ( ir >  +  ic > )  eq  x>igp-igp~l  ! same  value 
25  continue 

c**»*****************evaluate  f i 1  ter**************************** 

type  *,  '#  corrupted  pixels='.  ibp,  '  ^corrected  pixels='.  icp 
type  *.  '#  corrected  noise='.nbp.  '  ddisturbed  pixels='.igp 
so=0  O  ! Sum  of  orig 
sc=0.  0  ‘Sum  of  filtered 
poc=00  'cross  prod 
ssqo=0.  0  !  sum  orig  squared 
ssqc*0  0  !  sum  filt  squared 
dsum=0  0  !  sum  of  diff 
dsq=0  O  1  sum  dif  squared 
do  18  ir-31,482 
do  18  i c=31 i 482 

i p=dd < r < i r ) ♦  i  c )  'original  pixel 
if(ip  It  0)ip  =  ip+25<i 
po=f loa t ( ip+300 ) 

pc  =  f 1  oat ( cc <r (  i  r )  +  ic ) )  'filtered  pixel 

so=so+po  '  sum  of  orig 

sc=sc+pc  'sum  of  filtered 

poc =poc +po*p c  '  sum  cross  product 

ssqo=ssqo+po*p o  'sum  orig  squared 

s sqc =ssqc *pc *p c  'sum  filtered  squared 

dsum*dsum+(po-pc )  'sum  of  diffs 


18  dsq=dsq+<po-pc >*(p o-pc )  'sum  of  dif  squared 

c**************c omp u t e  &  display  results  of  eva 1 uat l on******************f « * 


dmean  =  d sum/g  'mean  of  difference  of  original  and  filtered  images 
dsd=sqr t < ( d sq/g ) -d mean*dmean >  ' stan  dev  of  difference 
corr«(g*poc-so*sc)/(sqrt< <g*ssqo-so*so)*(g*ssqc-sc*sc ) )  ) 
type  *.  'CORRELATION (orig inal. f 1 1 tered )= '.  coi  r 
type  *.  'Dif f (F ILTERED-OR IG >  Mean='. dmean 

type  *,  'Dif f (FILTERED-0R1G)  Si gma=  ' . dsd,  '  * 

type  *.  'Neighborhood  Replacement  FILTER  NOISF.=  '.p 
type  *.  'Gaussian  NOISE  Mean='.gm,  '  NOISE  Std  Dev='-sd 
call  exit 
end 

include  'mathfunc  for/nolist' 


program  onf4 

c  this  program  raads  daedalus  images  from  disk  then  introduces 

c  normal  noise  then  applies  criss-cross  filter  to  remove  the  noise 

c  then  computes  correlation  between  original  and  filtered  images 

include  'teamdef . f or/nol ist ' 

INTEGER *2  cc (262144)  ! corrupted  image 

integer*4  seed. r < 512) . rt (4) 

real*8  dsq,  dsum.  corr.  so.  sc.  poc.  ssqo.  ssqc 

byte  dd<262144)  (original  image 

character*20  f  (file  name 

data  rt/512.  1.  513,  511/ 

C************************* ********************************* 
seed-77737 
do  1  1  =  1,512 

1  r < i )=< i-1 >*512  (build  row  table 

C************************* ************************************** 

i=GT_INITIAL< 1,  2)  (init  term  function 

call  GT_PROMPT< '0NF4>  ')  (terminal  prompt 

c*************************************************************** 

88  i=gt_word( 'Disk  Fi 1 e  Name ( SFn ( DED) .  IMG) = f ( 1 : 20) > 

open (unit =7. fi le=f , for m= 'unformatted',  status='old', err =88) 
read (7, err=B8> , dd  (GET  IMAGE 
type  *.  '  IMAGE=  ',  f 
close! uni t=7 ) 

C ************************* ************************************** 
do  7  ir=23, 490  (converting  image  from  byte  to  i*2 
do  7  ic=23, 490 

ip=dd <r < ir )+ic >  (orig  byte 

if<ip  It  0)ip=ip+256  (  if  128-255  make  positive 
ip=ip+300  (shift  up  to  get  away  from  0  mean  gausian  noise 
7  cc (r ( ir >  +  ic )=ip  (put  in  contaminated  image. tho  not  yet  contam 

C ************* ************ **************************************** 
i=gt_real ( '#  Sigmas  for  bad  pisel  cr iter ion= ', sig ) 
i=g t_integ < 'Ne ighb orhood  Size  (ie  3x3=1;  5*5=2)=', is) 
i=gt_integ < '0=  Bad  Pit  Check;  1=  Skip=',ik) 
fn=( is+is*l )*< is+is+1 )  !n  for  neighborhood 
g=204304. O  1  #  of  pixels  in  452*452  image 
ibp=0  !#  pixels  contaminated 

if < g t_ask < 'Want  to  introduce  noi se  ', none >.  eq. no > goto  67 
i=gt^_real  <  'Enter  Percent  c on taminated  ' .  p  ) 
i=gt_real <  'Enter  Gaussian  Mean'.gm) 
i=gt_real ( 'Gaussian  std  dev.  ',sd> 

C ************* *****CONT AM I  NATE  WITH  NOISE************************ 

do  66,  ir=23,  490  (this  loop  enter  1-p  7.  gaussian  noise  mean  0,  sigma=sd 

do  66.  ic=23, 490  (corrupt  whole  picture 

if(ran(seed).  g  t  1  0-p)then 

cc(r(ir)  +  ic )=nint(gauss(gm, sd) )  (  noise 
if(ir  ge  31  and.  ir  le,  482)then 

if(ic  ge.  31  and  ic  le  482) ibp=ibp+l  !tt  pixels  corrupted 
end  i  f 
end  if 

66  continue 

C ****************************************************** ********* 

67  dsum=0  O 
dsq=0  0 

do  8  i r=31 .  482 
do  8  ic-31, 482 
ip=dd(r(ir)+ic ) 
if(ip  It  0)ip=ip+256 

pd  =  f 1  oat < i p+300 ) -f 1  oat < c c (r < ir > +  ic > >  (dif  orig-contam 
d sum=d sum+p d  1  Sum  of  diffs 


v.NV.y 


e 


dsq=dsq+pd*pd  ! sum  of  diff  squared 

dmean*dsum/g  (mean  of  difference  of  orig.  and  contaminated  images 
dsd^sqrt ( ( dsq/g )-dmean*dmean )  I  stan.  dev.  of  difference 
type  *, 'Diff (CONTAM. -ORIS  >Mean= dmean 
type  *.  'Diff (CONTAM. -ORIG.  >Sigma='.dsd 

C****************F I LTER************** ************* ***************** 
icp=0  ! #  pixels  corrected 
nbp=0  !#  noise  corrected 
igp»0  !#  good  disturbed 
ip  =  l 

9  do  25  ir=28+ip-l.  485-ip  +  l  (detect  bad  pixels  and  correct 

do  25  ic=28+ip-l. 485-ip+l 
x«float(cc(r<ir)+ic>> 
if ( i k  . ne.  0 ) goto  23 
dsq=0.  0 
dsum=0. 0 

do  21  iro--is.is  ! fn  window  compute  mean  and  sigma 
do  21  ico=-is< i s 

p  c  =  f  loat  (  cc  ( r  {  ir+iro  >■*  i  c +i  c  o  >  >  ‘get  pixel 
dsum-dsum+pc  (sum 

21  dsq=dsq+pc*p c  (sum  squared 

cmean-dsue/fn  (compute  local  mean 

csd=sqr t < dsq/fn-cmean*cmean )  (local  stan. dev 

*  f  <  <  x .  ge  cmean-si g*c  sd  >  and.  (x.  1  e  cmean+s i g*c  sd )>goto  25  (good 

23  p2=f loat ( cc ( r ( ir )+ic-r t ( ip > > > 

p3=float(cc(r(ir)+ic+rt(ip> )  ) 
x  x  =  x  ! save 

if(p2  .  le.  p3)then  (  x  p2  p3 

if(p2  . ge.  xlgoto  117  1  x  p2  p3 
if(x  .  le  p3)goto  25  ! p2  x  p3 

116  x=p3 

else  (x  p3  p  2 

if (p3  .  ge  xlgoto  116  !  x  p3  p2 
i f ( x  .  le.  p2>goto  25  (  p3  x  p2 

1 17  x=p2  (  p3  p2  x 
end  i  f 

icp=icp+l  (count  corrected  pixels  subtract  from  contaminated 
cc(r(ir)+ic)=nint(x> 
if  <  x  x  It.  100.  0 )  then 
nbp=nbp+l 
else 

if(x  ne  xx)igp=igp+l  ‘pixel  disturbed 
end  i  f 

25  continue 

type  *,  '#  corrupted  pixels='.ibp.'  #  corrected  pixels='iicp 
type  *,  '#  corrected  p i x e 1 s=  ' . nbp ,  '  #  disturbed  pixels*',  igp 
ip=ip+l 

if<ip  It.  5)goto  9 

C********************eval ua te  resu 1 ts ************* *************** 
so=0  0  (Sum  of  orig 
s c=0.  O  'Sum  of  filtered 
poc=0  0  (cross  prod 
ssqo=0  O  ( sum  orig  squared 
ssqc=0. 0  (sum  filt  squared 
dsum=0  0  ‘sum  of  diff 
dsq=0.  0  (sum  dif  squared 
do  18  i r=31 . 482 
do  18  i c  =3 1 . 482 

i p =dd < r < i r ) ♦ i c  )  ‘original  pixel 
if(ip  .  It  0>ip=ip+256 
po=float< l p+300 ) 
pc=float(cc(r(ir)+ic)) 
so=so*po  !  sum  of  orig 
sc=sc+pc  (sum  of  filtered 
poc“poc+po*pc  ‘sum  cross  product 
ssqo=ssqo+po*p o  ‘sum  orig  squared 
ssqc=ssqc+p c*p c  ‘sum  filtered  squared 
d su**d sum* ( p o-p c )  1  sum  of  diffs 

18  dsq=d sq-Mpo-pc  )*( p o-pc  )  (sum  of  dif  squared 

dmean=dsum/g  (mean  of  difference  of  original  and  filtered  images 

dsd=sqrt ( ( dsq/g )-dmean*dmean )  ‘stan  dev  of  difference 

corr*(g*poc-so*sc )/(sqrt( (g*ssqo-so*so)*(g*ssqc-sc*sc ) ) ) 

type  *.  ' CORRELATION  or  l  g  i  na  1 .  f  1 1  ter  ed  )  =  '.  c  orr 

type  *.  'Diff (FILTERED-ORIG)  Mean='.  dmean 

type  *,  'Diff (FILTERED-ORIG)  S i gma=  ' , d sd .  ' 

type  *, 'CRISS-CROSS  FILTER  NOISE*'. p 

type  *. 'NOISE  Mean='.gm,  '  NOISE  Std  dev='.sd 

call  exit 

end 

include  ‘mathfunc  for/nolist' 

C-9 


c  biweighted  3x3  filter  compute  weighted  mean  from  input  arvag 

subroutine  b  iuiei gh  t  (med.  xhat ) 
real*4  med(*).xhat 

c  H  of  elements  in  med=9, med ian-med ( S> 

sp=(med <7 > -med < 3> > / 1 . 349  !  compute  standard  deviation  from  quarti 
if < sp.  eq.  0) sp=l.  4B3  leant  have  0  sd 
s=*sp*5.  0 

c  calculate  initial  xhat  value 

sl=0.  0  I  sum 
•  2=0.  0  I  n 

do  lO  i  =  l<9  I  throw  out  outliers  >  5  std  devs 
if <abs (med ( i )-med <5) )  .  gt.  s)goto  10 
s l=s 1+med < i )  I  sum  of  non-outliers 
s2=s2+l.  0  In  of  non-outliers 

10  continue 

xhat=sl/s2  I init  computed  mean 
s=s+sp  I  <>  0  std  dev 

sl=0.  0  I  sum  of  weighted  non-outliers 
s2=0.  0  I  sum  of  weights 
do  30  i =1 . 9 

x=(med(i)-xhat)/s 
x  =  x*x 

if  (  x.  gt.  1.  0)  x  —  1  0 
wp  =  < 1 .  0- x  >»< 1  0- x ) 
s l=s 1 +wp*med ( i ) 

30  s2=s2«-wp 

xhat=sl/s2 

end 

********************************  ***■«■***«**#*■#  ****«■«■< 

REAL  FUNCTION  GAUSS (MEAN,  DEV) 

*  computes  a  random  number  from  a  normal  (gaussian)  distribution  with 

*  the  given  mean  and  deviation 

real  mean, dev, ran.  sum,  xx,  set_gauss_seed 
integer  Iran, i, idum 
data  iran/-l/ 


I 


«* 


c 

c 

c 


10 


20 


JO 


sum-0  0 
do  i " 1 . 50 

sum-tium+ran  (  i  ran  > 
enddo 

suro-sum/50  0 

xx=sqrt<12  0  •  t>0  0)  *  (sufn-0  5) 

gau5s-n*dev  mpan 

return 

entry  'jet_gaus sjped  (  j  d urn ) 
iran-i d  urn 

set  gau‘>‘,_5eed=0  0 

return 

end 

shell  sort 

sort  the  real  array  elements  *  <  l  )  to  *<n)  > n 

ascending  order 

subroutine  she  1 1  so r t ( x • n ) 
r ea 1 *4  x <  * ) ,  temp 
integer *4  i . n, nde 1 t j 
logical  inordr 

nd  v 1 ta~n 

if(ndelta  gt  1 ) t h e 1 1 
nde 1 ta=nde 1  to/2 
inordr-  true 
do  30/  j  r  1  /  n-ndpj  t»i 

i  f  (  x  (  i  )  g  t  x  <  i  -»  n  d  e  1 1  a  )  )  thou 
temp-  x  < l ) 
x  (  i  )  -  x  (  i  ♦  n  cl  v  1  t  a  ) 

*(  i  *-n delta)  te/np 
inordr2*  false 
end  i  f 
conti  nue 

if<  not  inordr  >  ooto  I'O 
goto  10 
end  i  f 
end 
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program  onfb 


c 

c 

c 

t 

include  'teamdef. for/nolist' 

INTEGER'S  cc (262144), sav, ino  'this  is  corrupted  image 
real  med<9) 
real*8  dsum,  dsq, 
integer*4  set*d,r(S12) 

byte  dd  (262144) ,  bb  (262144 )»  s  (2) .  s*>  (2)  !  original  image,  blank  image 
charac.ter*20  f  !  file  name 
equivalence  (s, sav), (inu, ss) 

c  *■**#■***•♦■»••*■*•**■**#*  * **#•*****«■  It '«#****»  V.  *K  e-K  *r,*-»***-iHHHHH*"»**«.*** 

seed=77737 
do  1  i=l, 512 

1  r  (  i  )  =  (  i-1  )*-512  .'build  row  table 

C*********-***##'**************-*********-#*  *•**■*  ****************** 
do  i  =  l, 262144  ! clear  noise  mask 
bb ( i )=0 
end  do 

i=GT_INITIAL( 1,  2)  fiviit  term  function 

call  GT_PROMPT < 'ONFB>  ')  !  terminal  prompt 

call  ipinit  ! initiate  deanza 

C  a-*-****-*-**-*-**##*****#-**#****#***  -tut  *■#  ns-  «.<•*  *•»*«•#*«■**#*•»***•»**  *  *  ** 

88  i=gt_word( 'Disk  File  Name  ( SFn  (DED).  IMG)*  ',  f  ( 1 :  20)  > 

open(unit=7, file=f , f orm= 'unf ormatted status='old err=8B) 

read (7, err=88> , dd  ‘GET  IMAGE 

close(unit=7) 

call  display (dd,  262144,  511, 0,  O)  {display  raw  image 
C********************************************************** 
do  7  ir=26, 487  !  converting  image  from  byte  to  i*2 
do  7  ic=26, 487 

ip=dd(r(ir)+ic )  ! orig  byte 

if (ip  .It.  O) i p=i p+256  !  i f  128-255  make  positive 
ip=ip+300  {shift  up  to  get  away  from  O  mean  gausian  noise 
7  cc (r ( ir )+i c >=i p  {put  in  contaminated  image, tho  not  yet  contam 

C***#***»*##**»*»**»»#**** *»***«« ********  **********  ********** 
i=gt_real('#  Sigmas  for  bad  pixel  cr  iterion® ', sig > 
i=gt_integ ( 'Neighborhood  Size  (ie  3x3=1 »  5x 5=2 )=',  i  s  ) 
fn= ( is+is*l )*(is+i s+1 )  !  n  for  neighborhood 
g =204304. 0  ! #  of  pixels  in  452*452  image 
ibp=Q  !tt  pixels  corrupted 

if (gt_ask ( 'Want  to  introduce  noi se ', none ) .  eg.  no ) goto  67 
i=gt_integ ( 'Color  of  Noise®', ino) 

**9^_ra«‘l  1  'Enter  Percent  contaminated  ( 10‘/.=  .  l)',p) 
i=gt_real ( 'Gaussian  mean=',qm)  {allow  any  gaussian  mean 
i=gt_real  (  'Gaussian  std.duv.  '.sd) 

c*********************Contaminate  with  noi so**********-************* 
sav=254 

do  66,  ir=26, 487  {this  loop  enter  1-p  '/.  gaussian  noise  mean  0,  sigma- sd 

do  66,  ic=26, 407 

if(ran(seed)  .  gt.  1.0-p)then 

cc  (r(ir)-*-ic)=nint(gatiss(gm,  sd)  >  !  noise 

dd ( r (  ir )  +  i c )  =  ss ( 1 )  {put  noise  in  picture  color  selected  (ino) 
bb (r ( ir )+i c )=s ( 1 )  ! put  noise  in  blank  image  (mask)  (ino) 
if(ir  .  ge.  31  .and.  ir  .  le.  482)then 

if(ic  .ge.  31  .and.  i <•  .  le.  402 )ibp  =  ibp  +  l  !#  pixels  corrupted 
end  i  f 
end  if 
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66  continue 

C**********************d  i  f  between  or i g/ c  ontaminat ed******************** 

67  icp=0  !#  pixels  corrected 
igp=0  !#  good  pixels  fixed 
nbp-0  I H  bad  pixels  fixed 

call  d i sp lag < dd. 262144,  31 1 , 0/ 0>  ! display  noise 

if  <g  t_ask  <  'Disp  lag  Blank  '*  none),  eq.  YE3)call  disp  lay  (bb.  262144.  81 1*  0<  0) 
c****#**#*#*****FILTER  IMAGE  a##***#********#****##****###*###***## 
do  25  ir=31.482  ! detect  bad  pixels  and  correct 
do  25  ic=31.  482 

x^floaticcfriir-J+ic)  >' 
dsq=0.  0 
dsum=0. 0 

do  21  iro=-iSiis  !  11x11  window  compote  mean  and  sigma  for  windoui 
do  21  ico=-is,  is 

pc=f  loat  (cc  <  r  (  ir+ir  i>  )  +  ir.  +  ico  )  )  (get  pixel 
dsum=d sum  i-g  e:  !  sum 

21  d sq=d sq+pc *p c  !  sum  squared 
cmean--d num/f r>  (compote  local  mean 
csd~sqrt < dsq/fn”cmean*cmean )  (local  s tan. dev 

if((x.  go.  cmean-si  g*c  sd  ).and.  (  x.  le.  cmean+s  i  g*c  s  d  >  )  goto  25  (  good 

icp=icp+l  (count  corrected  pixels 

ict=0 

do  22  iro--l.l  (extract  3x3  windoui 
do  22  ico=-l.  1 

pc=f loat (cc (r ( ir+ivo )+i c+ico ) )  (get  pixel 
ic t«ict+i 

22  med(irt)=pc 

call  shellsortimedj 9> 

sav=ni nt (med ( 5) ) -300  (median  pixel 

if(sav  .  It.  0)sava0 

dd  (r  ( ir  )  +  ic  >=«;<  1 )  (corrected  pixel  in  orig  pic 
i f < cc ( r ( ir )+ ic )  .It.  20O)then  (bad  pixel  fixed 
nbp=nbp+l  (bad  pixel  fixed 
bb  (r  <  ir  )  +  i  c  )-<>  (black  pixel 
else 

igp-:igp  +  l  (good  pixel  messed  op 
bb (r (ir )+ic  )=s( 1 )  (yellow  pixel 
end  i  f 

25  continue 

c  ft****************-* -.H.-eval  uate  f  i  1  ter*********#- ****************** 
call  d i sp lay ( d d. 262144. 51 1 . 0. 0 )  (display  corrected  pix 
type  *.  '#  corrupted  pixels^'.  ihp.  '  41c orrec ted  pixelsB'i  icp 
type  *»  '#  corrupted  pixels  fixed”'*  nbp»  '#guuil  made  bcui-'.igp 
i f  ( g  t_ack  ( 'Di  sp  lay  Bl. ink '»  none: ).  eq.  YES )  then 
call  displayibb.  262144.  511.0.0' 
end  i  f 

c**************c omp ute  &  display  results  of  evaluation********************* 
type  *,  'FILTER  NOISE*'/ p 

type  *.  'Gaussian  NOISE  Mean”' >  gin*  '  NOISE  Std  Dev-  sd 

if <gt_ask ( 'Quit none),  eq.  YESlcall  exit  (gives  time  to  display  image 
end 

include  'mathfunc.  for/nol i  u  l  ' 
include  'shellsort.  for/nul ist ' 
include  'ipinit. for/nolist ' 


